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INTRODUCTION 


1.  GENERAL 

Performance  prediction  for  a  ramjet  propulsion  system  requires 
a  detailed  analysis  of  the  flow  through  each  of  its  major  components: 
the  inlet,  the  combustor,  and  the  exhaust  nozzle.  The  objective  of 
this  research  was  to  develop  a  procedure  for  calculating  the  flow 
field  within  the  converging-diverging  nozzle  of  a  high-speed  (Ms6) 
ramjet  propulsion  system. 

Initial  data  (at  the  nozzle  inlet)  for  the  nozzle  flow  field 
analysis  is  provided  by  the  results  of  the  combustor  analysis.  Ramjet 
combustors  are  generally  axisymmetric  in  shape  with  slowly-varying  or 
constant  cross-sectional  area.  Techniques  for  analysis  and  performance 
prediction  for  a  high-speed  ramjet  combustor  are  still  under  develop¬ 
ment.  One  technique  [Reference  ( 1  )]  is  based  on  a  one -dimensional 
compound  flow  analysis  of  the  combustor  flow  field.  The  flow  is  broken 
into  a  number  of  coaxial  flow  streams  and  the  one-dimensional  flow 
equations  are  solved  within  each  stream.  The  static  pressure  is 
matched  radially  across  each  of  the  streams  resulting  in  a  uniform 
static  pressure  distribution  across  the  combustor  exit.  The  fuel 
injection  processes,  fuel  droplet  dynamics,  atomization,  mixing,  vapor¬ 
ization,  and  combustion  kinetics  are  all  included  in  the  flow  model. 
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The  final  result  of  the  combustor  analysis  is  a  complete  set  of  known 
flow  properties  and  chemical  species  concentrations  at  the  nozzle 
entrance  where  the  converging  geometry  begins  to  introduce  two- 
dimensional  flow  effects.  The  flow  at  the  combustor  exit  will  not 
necessarily  be  in  chemical  equilibrium.  At  high  speeds,  the  stagna¬ 
tion  temperature  within  the  flow  will  be  large  and  the  effects  of 
finite-rate  chemical  kinetics  become  important.  Also  residence  times 
within  the  combustor  may  not  be  sufficient  for  the  flow  to  reach  chem¬ 
ical  equilibrium.  The  technique  developed  in  the  subject  research 
will  accept  a  nonuniform,  nonequilibrium,  distribution  of  initial  con¬ 
ditions  at  the  nozzle  inlet. 

The  analysis  within  the  nozzle  can  be  divided  into  two  parts: 
the  analysis  of  the  subsonic  and  transonic  flow  fields  in  the  axisym- 
metric  nozzle  entrance  and  throat  regions,  and  the  analysis  of  the 
supersonic  flow  field  in  the  nozzle  divergence.  The  effects  of  finite- 
rate  kinetics  must  be  included  in  both  parts  of  the  analysis.  If  the 
nozzle  divergence  is  axisymmetric.  Reference  (2)  can  be  used  for 
analysis  of  the  supersonic  flow  field.  If  the  divergence  is  not  axi- 
symnetric,  due  to  integration  with  the  airframe,  then  a  three- 
dimensional  analysis  is  required.  A  highly  accurate,  three-dimensional, 
method-of -characteristics  scheme  for  supersonic  flow,  including  finite- 
rate  chemical  kinetics  effects,  has  been  developed  by  Cline  and  Hoffman 
(3).  The  subject  research  has  been  concerned  with  the  analysis  of  the 
subsonic  and  transonic  flow  regions.  Both  axisymmetric  and  planar 
nozzle  geometries  are  considered  and  the  nozzle  may  have  a  centerbody. 
Though  the  present,  discussion  focuses  on  ramjet  propulsion,  the 
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technique  which  is  described  is  applicable  to  any  subsonic-transonic 
flow  field  analysis  where  the  effects  of  finite-rate  chemical  kinetics 
must  be  considered.  The  computer  program  which  evolved  from  this  work 
is  an  extension  of  that  developed  by  Cline  (4). 

2.  LITERATURE  REVIEW 

The  literature  surveyed  for  this  research  can  be  divided  into 
three  groups:  1)  one-,  two-,  and  three-dimensional,  steady,  chemically 
reacting  nonequilibrium  flow,  2)  two-dimensional,  unsteady,  nonreacting 
flow,  and  3)  unsteady,  chemically  reacting,  nonequilibrium  flow.  Each 
of  these  groups  is  considered  in  turn  in  the  following  discussion. 


Steady,  Chemical ly  Reacting,  Nonequilibrium  Flows 

There  is  a  very  large  amount  of  literature  concerning  this  topic. 
Most  of  the  literature  deals  with  one-  and  two-dimensional  flows.  For 
the  one-dimensional  case,  the  solutions  can  be  categorized  as  analyti¬ 
cal,  semiempirical ,  and  numerical.  Analytical  solutions  are  possible 
only  for  the  very  simplest  nonequilibrium  processes  and  often  approxi¬ 
mations  such  as  linearization  and  the  assumption  of  a  high  Mach  number 
are  introduced.  Reference  (  5  )  provides  examples  of  analytical  solu¬ 
tions.  Perhaps  the  most  often  used  semiempirical  method  is  the  "sudden 
freezing"  approximation  (  5)  where  equilibrium  flow  is  assumed  upstream 
of  an  empirically  chosen  "freezing  point,"  and  frozen  flow  is  assumed 
downstream  of  that  point.  Accurate  numerical  solutions  are  also  avail¬ 
able  and  have  been  computed  by  a  number  of  investigators.  Examples 
include  the  work  by  Emanual  and  Vincenti  (6  ,  7  ),  Sarli  (8)»  Lordi 


et  al .  (9),  and  Frey  et  al.  (10).  Reference  (10)  is  widely  used  by  the 
propulsion  industry. 

Two-dimensional,  steady,  supersonic  flows  with  chemical  nonequili¬ 
brium  have  been  computed  by  Widawsky  ( 1 1 ) ,  Zupnik  et  al.  ( 12),  Craig 
(13),  Bartlma  (14),  and  Kliegel  et  al.  (2).  Reference  (2)  considers 
the  elements  C,  H,  0,  N,  F,  and  Cl,  19  gaseous  species  containing 
those  elements,  and  48  chemical  reactions  that  may  occur  between  those 
species.  A  method-of -characteristics  procedure  is  used.  The  super¬ 
sonic  initial-data  line  required  to  start  the  characteristic  calcula¬ 
tions  can  be  input  or  it  can  be  calculated  by  a  transonic  analysis 
within  the  program.  The  characteristic  equations  governing  the  fluid 
dynamic  variables  are  integrated  by  a  second-order  (modified  Euler) 
explicit  method.  The  species  continuity  equations  are  integrated  using 
a  first-order,  implicit  method.  It  is  significant  to  note  that  the 
authors  of  Reference  (2)  state  that  even  small  interpolation  errors 
in  species  concentrations  cause  severe  stability  and  accuracy  problems 
in  the  numerical  integration  of  the  species  continuity  equations.  This 
fact  will  be  discussed  more  fully  in  Section  III. 

Cline  and  Hoffman  (3)  extended  an  isentropic  flow  three- 
dimensional  method-of -characteristics  scheme,  developed  by  Ransom  et 
al.  (15),  to  include  the  effects  of  finite-rate  chemical  kinetics. 

Their  analysis  uses  the  same  chemical  kinetics  model  as  used  in  Refer¬ 
ence  (2).  Initial-value  plane  data  can  be  obtained  from  the  results 
of  an  analysis  using  Reference  (2),  but  in  that  case,  no  account  is 
made  for  two-dimensional  effects  upstream  of  the  throat.  The  subject 
research  provides  a  two-dimensional,  nonequilibrium  initial-value 


plane  for  use  in  the  analysis  of  Reference  (3).  Cline  and  Hoffman  (3) 
note  the  importance  of  using  correct  initial  data.  They  found  that 
incorrect  data  may  result  in  thermal  compressions  and  shocks  due  to  the 
rapid  change  in  chemical  composition  just  off  the  initial -value 
surface.  As  in  Reference  (2),  the  fluid  dynamic  equations  are  inte¬ 
grated  using  a  second-order  explicit  method,  but  the  species  continuity 
equations  are  integrated  using  a  second-order,  implicit  scheme.  This 
particular  scheme  was  chosen  after  an  extensive  review  of  various 
explicit  and  implicit  methods  for  integrating  the  species  continuity 
equations. 


Unsteady,  Nonreacting  Flow 

The  time-dependent  technique  for  the  solution  of  steady-state 
converging-diverging  nozzle  flows  has  been  used  by  a  number  of  inves¬ 
tigators  [References  (4)  and  (16)  to  (23)].  References  (16)  to  (22) 
provide  generally  good  results  for  typical  problems  but  the  techniques 
yield  relatively  long  computational  times.  Cline  (23)  attributes  this 
to  "inefficient  numerical  schemes  and  poor  treatment  of  the  boundaries 
resulting  in  the  requirement  for  excessively  fine  computational  meshes." 
Cline  (4,  23)  uses  MacCormack's  method  for  interior  points  with  the 
governing  equations  of  motion  in  nonconservation  form.  Implicit  arti¬ 
ficial  viscosity  is  present  within  this  method  so  explicit  artificial 
viscosity  is  added  only  for  flows  with  shock  waves.  The  boundary 
points  are  computed  using  a  reference-plane  characteristic  scheme. 
Solutions  to  representative,  inviscid,  nozzle  problems  are  obtained 
in  approximately  one  minute  using  a  CDC  6600  computer. 


Vamos  and  Anderson  (24)  have  used  a  time-dependent  scheme  to  com¬ 
pute  the  one-dimensional,  nonequilibrium  flow  of  reacting  gas  mixtures. 
They  cite  several  advantages  of  their  approach  as  compared  to  steady, 
one-dimensional  analyses: 

1)  the  use  of  relatively  few  grid  points  with  relatively  large 
spatial  increments  through  the  nozzle,  including  the  near 
equilibrium  subsonic  portion  of  the  nozzle. 

2)  the  straightforward  manner  in  which  completely  nonequilibrium 
subsonic,  supersonic  flow  can  be  treated.  In  particular,  the 
usual  singularity  in  the  throat  region  associated  with  steady 
flow  analyses  is  not  present  in  the  unsteady  approach. 

MacCormack's  method  (which  is  explicit)  was  used  in  this  work  for  all 
the  governing  equations  (Including  the  species  continuity  equations). 
Very  long  computation  times  were  reported  (M  hour  on  a  CDC  6400)  but 
the  results  were  quite  good.  No  general  computer  programs  have  been 
developed  to  treat  multidimensional,  reacting  flows  by  an  unsteady 
approach. 

3.  REMARKS  ON  COMPUTATIONAL  TIME 

Stability,  accuracy,  and  computational  time  are  generally  the 
most  Important  factors  used  In  determining  the  merit  of  a  particular 
numerical  method.  In  the  subject  research,  computational  time  has 
been  especially  Important  for  several  reasons.  First,  the  flow  field 
to  be  modeled  in  the  nozzle  entrance  and  throat  regions  is  a  two- 
dimensional,  subsonic  and  transonic  flow.  For  steady  subsonic  flow. 


the  governing  partial  differential  equations  are  elliptic  while,  for 
steady  supersonic  flow,  those  equations  are  hyperbolic.  In  the  tran¬ 
sonic  flow  regime,  the  equations  change  from  elliptic  to  hyperbolic 
(i.e.,  they  are  mixed).  Techniques  for  solving  mixed  flow  problems 
are  not  well  developed,  but  hyperbolic  problems  have  been  studied  and 
solved  for  many  years.  The  mathematical  difficulties  associated  with 
a  mixed  flow  problem  can  be  avoided  by  solving  the  unsteady  flow  equa¬ 
tions  which  are  hyperbolic  in  both  subsonic  and  supersonic  regions. 

Then  the  steady  flow  solution  can  be  obtained  as  the  asymptotic 
solution  to  the  unsteady  equations,  with  steady  flow  boundary  conditions 
applied,  for  large  time.  The  disadvantage  of  this  approach  is  that 
another  independent  variable  is  introduced  into  the  analysis  and  the 
unsteady  solution  may  have  to  be  advanced  through  many  time  steps  to 
achieve  the  steady  state.  If  the  allowable  time  step  is  small,  and/or 
the  time  to  compute  each  solution  plane  is  large,  then  long  computa¬ 
tional  times  can  result.  The  time  step  size  is  determined  from 
stability  considerations  for  the  given  equations  and  numerical  methods 
used.  Introducing  finite-rate  chemical  kinetics  into  the  analysis 
significantly  increases  the  time  required  to  compute  each  solution 
plane.  This  is  because  the  number  of  equations  in  the  mathematical 
model  increases  by  the  number  of  chemical  species  considered,  and  the 
calculations  required  to  compute  the  species  source  function  and  to 
integrate  the  species  continuity  equations  are  very  time  consuming. 

The  nature  of  the  species  continuity  equations  must  also  be  con¬ 
sidered  in  a  discussion  of  computational  time.  For  flows  near 


equilibrium,  the  species  continuity  equations  are  "stiff,"  and  as  a 
result,  they  are  difficult  to  solve  numerically  (see  Appendix  D). 
Standard  explicit  integration  techniques,  when  applied  to  stiff 
differential  equations,  are  unstable  except  for  very  small  step  sizes. 
This  problem  occurred  in  the  work  of  Vamos  and  Anderson  (24).  The  use 
of  implicit  methods,  however,  removes  the  stability  problems  and  allows 
more  reasonable  time  steps.  In  addition  to  obtaining  a  solution  to 
the  subject  problem,  a  significant  objective  of  this  research  has  been 
to  achieve  the  solution  in  reasonable  computational  times. 


SECTION  II 


MODELING  THE  PROBLEM 


1.  PHYSICAL  MODEL 

In  Section  I  it  was  noted  that  the  initial  data  for  the  nozzle 
flow  field  analysis  is  provided  by  the  results  of  the  combustor 
analysis.  Within  the  combustor,  the  fuel  injection  processes,  fuel 
droplet  dynamics,  atomization,  mixing,  vaporization,  and  combustion 
kinetics  may  all  be  considered.  However,  by  the  time  the  flow  has 
reached  the  combustor  exit,  it  is  assumed  in  this  analysis  that  it  is 
a  mixture  of  thermally  perfect  gases  with  no  condensed  phases.  The 
distribution  of  fluid  dynamic  properties  and  the  species  concentrations 
may  or  may  not  be  uniform  across  the  exit.  Also,  the  flow  may  be  in 
chemical  nonequilibrium  at  the  combustor  exit.  The  flow  within  the 
nozzle  is  treated  as  continuous,  inviscid,  and  adiabatic.  Body  forces 
are  neglected.  No  provisions  have  been  made  for  diffusion  or  turbu¬ 
lence  modeling  within  the  flow.  Finally,  only  chemical  nonequilibrium 
is  treated  in  this  analysis;  the  flow  is  assumed  to  be  everywhere  in 
instantaneous  translational,  rotational,  and  vibrational  equilibrium. 

2.  MATHEMATICAL  MODEL 

The  equations  which  correspond  to  the  physical  model  described 
above  are  given  in  Reference  (25)  and  are  discussed  in  Appendix  A. 
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These  equations  are  the  global  continuity  equation,  the  momentum 
equation,  the  energy  equation,  the  species  continuity  equations,  and 
the  thermal  and  caloric  equations  of  state.  In  vector  form,  they  are: 
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where  p  is  the  fluid  density,  V  is  the  velocity,  P  is  the  pressure, 
h  is  the  system  enthalpy,  T  is  the  temperature,  are  the  species 
densities,  are  the  species  mass  fractions,  are  the  species  source 
functions,  h^  are  the  species  enthalpies,  are  the  species  gas  con¬ 
stants,  Cp.|  are  the  species  specific  heats  at  constant  pressure,  and 
hj°  are  the  species  energies  of  formation. 

In  Appendix  A,  equations  (3)  and  (4)  are  manipulated  to  forms 
which  are  more  convenient  for  analysis  and  the  entire  set  of  governing 
equations  is  written  for  two-dimensional  axisymmetric  or  planar  flow. 


The  resulting  equations  are: 


Pt  +  UPx  +  VPy  +  pux  +  pvy  +  £(pv/y)  =  0  (?) 

ut  +  uux  +  vuy  +  Px/P  =  0  (8) 


(9) 

(10) 

(ID 


The  subscripts  in  these  equations  denote  partial  differentiation  and 
e  in  equation  (7)  is  zero  for  planar  flow  and  one  for  axisymmetric 
flow.  The  terms  a^  and  are  the  frozen  speed  of  sound  and  the  ratio 
of  frozen  specific  heats  respectively.  Note  that  the  species  conti¬ 
nuity  equations  (11 )  are  coupled  to  equations  (7)  to  (10)  by  the  energy 
equation  source  term  <1^. 

Since  the  interior  points  in  this  analysis  are  to  be  treated  by 
a  fixed  grid  technique,  it  is  convenient  to  transform  the  physical 
(x,y,t)  plane  to  a  rectangular  (C,n,x)  plane  in  which  the  differencing 
is  performed.  The  following  coordinate  transformation  is  used  (see 
Figure  1): 
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where  yc(x)  and  yw(x)  represent  the  nozzle  centerbody  and  wall  coor¬ 
dinates,  respectively,  as  functions  of  x.  Applying  this  transformation 
to  equations  (7)  to  (11)  yields  the  forms  of  the  continuity,  momentum, 
energy,  and  species  continuity  equations  which  are  used  in  this  analy¬ 
sis. 
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Note  that  equation  (17)  may  be  written  as 


(i=l,...,n) 


where  DC^/Dt  is  the  change  in  C..  following  a  particle  path  in  the  com¬ 
putational  plane. 


3.  CHEMICAL  KINETICS 

The  species  source  function  appears  in  the  forcing  terms  of  both 
the  energy  equation  (16)  and  the  species  continuity  equations  (18). 
Before  the  species  source  function  can  be  computed,  a  reaction  mechan¬ 
ism  must  be  specified.  The  reaction  mechanism  used  in  this  research 
is  the  same  as  that  used  in  Reference  (3)  with  the  addition  of  the 
hydrocarbon  reaction  discussed  in  Appendix  B.  Unburned  hydrocarbons 
may  be  present  in  the  gas  mixture  at  the  combustor  exit  and  therefore, 
provision  for  these  species  must  be  made  in  the  chemical  kinetics  model 
Edelman,  et  al.  (26)  and  Edelman  and  Harsha  (27)  have  developed  a 
kinetics  model  for  use  in  combustor  analyses  which  includes  a  "sub- 
global"  partial  oxidation  step: 


C  H  +  £  09  +  ?  H,  +  nCO 
n  m  c  c  c  t. 


Note  that  this  reaction  proceeds  in  the  forward  direction  only.  The 
reaction  rate  coefficient  for  equation  (19)  is  determined  empirically 
and  is  given  by  (Ref.  27): 


a.  =  6.0  x  1014  P0,3  c i  ..  Cn 


(20) 


where  P  has  units  of  atmospheres  and  c  represents  species  concentra¬ 
tions  in  g-moles/cc.  For  a  nonequilibrium,  chemically  reacting  flow 
of  a  system  of  thermally  perfect  gases  composed  of  the  six  elements 
carbon,  hydrogen,  oxygen,  nitrogen,  flourine,  and  chlorine.  Gold  and 
Weekly  (28)  have  shown  that  19  species  and  48  chemical  reactions  should 
be  considered.  Tables  B-I  and  B-II  of  Appendix  B  show  these  species 
and  reactions.  Experience  with  this  and  other  reaction  mechanisms  has 
indicated  that  there  may  be  cases  when  not  all  of  the  reactions  corres¬ 
ponding  to  a  given  set  of  chemical  species  should  be  included  in  the 
analysis.  The  number  of  reactions  in  the  reaction  mechanism  is  a  sig¬ 
nificant  parameter  affecting  execution  time.  Therefore,  no  more  reac¬ 
tions  than  are  necessary  should  be  included  in  the  reaction  mechanism. 

A  general  reaction  equation  which  is  used  to  represent  any  reac¬ 
tion  mechanism  is 
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where  v! .  and  v!.' .  are  the  stoichiometric  coefficients  of  the  reactants 

•  J  I  J 

and  products  respectively,  denotes  the  ith  chemical  species,  and  k^j 
and  k  .  are  the  forward  and  reverse  reaction  rate  coefficients,  respec- 

•  J 

tively,  for  the  jth  reaction  of  the  m  reactions  in  the  reaction 
mechanism. 

Reaction  rate  coefficients  are  not  readily  predicted  at  the 


present  time  and  therefore,  they  are  generally  found  experimentally. 
The  form  for  the  reverse  reaction  rate  coefficient  which  is  used  in 
this  analysis  is 


where  a..,  n.,  and  b.  are  empirical  coefficients,  R  is  the  universal 
1 J  J  J 

gas  constant,  and  T  is  the  local  gas  temperature. 

Expressions  for  the  species  source  function  for  both  dissociation 
recombination  reactions  and  binary  exchange  reactions  are  developed  in 
Appendix  B.  For  dissociation-recombination  reactions,  the  species 
source  function  is: 
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For  binary  exchange  reactions  the  species  source  function  is: 


In  equations  (23)  and  (24),  K0  *  is  the  equilibrium  constant  based  on 
partial  pressures,  Mj  is  the  third  body  reaction  rate  ratio,  are 
the  species  mass  fractions,  and  in.  are  the  species  molecular  weights. 

The  effect  of  different  third  bodies  on  dissociation-recombination 
reactions  is  accounted  for  through  the  use  of  third  body  reaction  rate 
ratios  M.  as  in  equation  (23).  Appendix  B  contains  a  discussion  of 

J 

these  ratios. 

Note  that  the  difficulties  in  making  accurate  measurements  for 
the  determination  of  reactions  rate  coefficients,  the  necessity  of 
extrapolating  experimental  data  from  one  situation  to  another,  and  the 
uncertainty  as  to  the  reaction  mechanism,  are  all  significant  factors 
limiting  the  accuracy  of  the  analysis  presented  in  this  research.  How¬ 
ever,  experience  with  one-dimensional  analyses  of  reacting  flows 
indicates  that  the  chemical  kinetics  model  used  here  is  adequate  for 
performance  prediction. 


SECTION  III 
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THE  UNSTEADY  AND  INCONSISTENT  SCHEMES 

1.  GENERAL 

The  governing  equations  presented  in  Section  II  can  be  thought 
of  as  being  composed  of  two  groups:  1)  the  equations  governing  the 
fluid  dynamic  variables  (i.e.,  the  global  continuity,  momentum,  and 
energy  equations),  and  2)  the  species  continuity  equations.  For 
convenience,  the  first  group  will  hereafter  be  called  the  fluid 
dynamic  equations.  Two  approaches  to  the  solution  of  the  complete 
set  of  governing  equations  have  been  considered  and  are  discussed 
below.  The  first  approach  is  called  the  "unsteady  scheme"  and,  though 
it  has  not  been  implemented,  it  is  helpful  in  understanding  the  second 
approach.  The  second  approach  is  called  the  "inconsistent  scheme"  and 
it  has  been  used  successfully  in  the  overall  numerical  algorithm  of 
this  research. 

2.  THE  UNSTEADY  SCHEME 

Consider  an  initial-data  surface  for  the  subject  problem  at  time 
level  N.  Figure  2  illustrates  the  unsteacjy  scheme  at  an  arbitrary 
grid  point  within  the  computational  mesh.  The  fluid  dynamic  variables 
uN,  v\  PN,  and  pN  and  the  species  mass  fractions  C^  are  known  at  all 
mesh  points.  The  superscripts  denote  the  time  level.  In  the  unsteady 
scheme,  the  solutions  of  the  fluid  dynamic  equations  are  advanced  one 
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unsteady  scheme. 


time  step  by  the  methods  of  Reference  (38)  to  give  uN+1,  v^+*,  P^+*, 
and  p^+*.  The  are  held  constant  during  this  step.  In  words,  the 
species  continuity  equations  given  by  equation  (18)  require  that  the 
total  change  in  mass  fraction  along  a  particle  path,  for  each  of  the 
chemical  species,  is  equal  to  c^/p.  Therefore,  using  the  solution  for 
uN+1  and  vN+1  and  the  time  step  size,  the  particle  path  is  traced  back 
from  the  grid  point  at  time  N+l  to  its  origin  in  the  initial  data  sur¬ 
face  (time  level  N).  Then,  by  interpolation  in  the  two-dimensional 

N  N  N  N  N 

initial-value  surface,  values  of  u*  ,  v*  ,  P*  ,  p*  ,  and  at  the 
origin  are  determined.  The  subscript  *  indicates  values  at  the  parti¬ 
cle  path  origin  which  will  generally  not  be  at  a  computational  mesh 
point.  Given  u*N,  v*N,  T*N,  p*N,  C.*N,  uN+1,  vN+1,  TN+1,  and  pN+1,  it  is 
possible  to  integrate  the  species  continuity  equations  forward  along 
the  particle  path  and  to  compute  C.^+*.  When  the  steady  state  is 
achieved,  the  same  particle  path  and  *  values  will  be  computed  for  each 
additional  time  step  and  therefore,  the  solution  will  not  change  with 
time.  There  are  a  variety  of  unit  processes  which  might  be  imple¬ 
mented  with  this  approach.  For  example,  one  might  predict  the  solution 
for  uN+1,  vN+1,  PN+1,  and  pN+1,  generate  the  particle  path  and  inter¬ 
polate  in  the  initial-data  plane,  integrate  the  species  continuity 

N-t-l 

equations  forward  to  give  C?  ,  and  then  repeat  these  same  steps  in 
a  corrector  procedure. 

The  unsteady  approach  suffers  several  serious  drawbacks.  Time 
consuming  interpolation  for  4+n  variables  in  the  two-dimensional 
initial-data  surface  is  required.  Also,  as  noted  in  Reference  (2), 
interpolation  for  species  concentrations  is  likely  to  produce  numerical 
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difficulty.  Values  of  the  forcing  term  in  the  energy  equation  for 
high  temperature,  near  equilibrium  conditions  are  found  to  be  very 
sensitive  to  variations  in  the  species  mass  fractions.  Therefore,  any 
interpolation  error  in  the  values  of  C.  could  produce  dramatic  errors 
i n  *  ^eca^  that  coup'es  the  fluid  dynamic  equations  and  the 
species  continuity  equations.  The  unsteady  approach  is  useful  because 
it  is  conceptually  straightforward  and  provides  a  clear  picture  of 
convergence  at  the  steady  state.  Also,  it  provides  a  basis  from  which 
the  inconsistent  approach  can  be  understood. 

3.  THE  INCONSISTENT  SCHEME 

Consider  the  ordinary  differential  equation 

$  85  f(y.t)  (25) 

The  finite  difference  form  of  this  differential  equation  is  "consistent" 
if  the  finite  difference  equation  approaches  the  differential  equation 
as  At  -*■  0.  The  overall  algorithm  used  in  this  analysis  is  inconsis¬ 
tent  because  of  the  treatment  of  the  species  continuity  equations. 
Specifically,  in  the  inconsistent  approach,  the  solutions  for  the  fluid 
dynamic  variables  are  advanced  one  time  step  from  an  initial  data  sur¬ 
face  in  the  same  manner  as  for  the  unsteady  approach  while  the  species 
mass  fractions  are  held  constant.  Then,  the  flow  at  time  level  N+l 
is  assumed  to  be  steady  and  as  many  as  four  streamlines  are  traced 
through  the  flow  field  with  origins  at  selected  locations  along  the 
nozzle  inlet.  The  species  continuity  equations  are  then  integrated 
along  these  streamlines  from  the  nozzle  inlet  to  the  exit.  This 
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process  of  advancing  the  fluid  dynamic  variables  and  then  integrating 
the  species  continuity  equations  along  streamlines  is  repeated  until 
convergence  at  the  steady  state  is  achieved.  This  approach  is  incon¬ 
sistent  because  the  finite  difference  form  of  the  species  continuity 
equations  does  not  approach  equation  (18)  as  At-*  0  (SC^/at  is 
neglected) . 

Figure  3  is  helpful  in  understanding  the  inconsistent  approach 
and  its  relationship  to  the  unsteady  approach.  Note  that  it  corres¬ 
ponds  to  the  special  case  of  v  equal  to  zero  and  a  streamline  passing 
through  computational  mesh  points.  Advancing  the  fluid  dynamic  var¬ 
iables  one  time  step  gives  values  of  u^+1,  v^+1,  P^+1,  and  pN+1  at 
point  C  in  Figure  3.  It  would  be  possible  to  trace  a  particle  path 
back  to  point  A,  corresponding  to  l-l  in  the  computational  mesh, 
where  values  of  u,  v,  P,  p,  and  Ci  could  be  determined  by  interpola¬ 
tion  (between  values  at  L-l,  M,  N  and  those  at  L-l,  M,  N-l).  Then, 
just  as  in  the  unsteady  approach,  the  species  continuity  equations 
could  be  integrated  forward;  in  this  case,  from  point  A  to  point  C. 

The  proper  time  step  for  integration  of  the  species  continuity  equa¬ 
tions  would  be  found  from  the  velocities  at  A  and  C  and  the  distance 
between  those  two  points.  This  time  step  would,  in  general,  not  be 
equal  to  the  time  increment  between  N  and  N+l.  In  the  inconsistent 
scheme,  the  species  continuity  equations  are  integrated  between  points 
B  and  C  with  the  time  step  determined  by  the  velocities  at  B  and  C 
and  the  distance  between  those  two  points.  If  the  values  of  u,  v,  P, 
p,  and  are  the  same  at  points  A  and  B,  then  the  inconsistent  approach 
gives  the  same  result  as  the  unsteady  scheme.  Also,  in  the  limit  of 
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Figure  3.  The  inconsistent  scheme. 


steady  state,  the  values  of  the  fluid  dynamic  variables  and  species 
concentrations  do  not  change  between  points  A  and  B. 

The  benefit  of  using  the  inconsistent  scheme  is  that  integration 
of  the  species  continuity  equations  always  proceeds  from  the  computed 
values  at  B  to  yield  new  values  of  the  species  mass  fractions  at  point 
C.  Interpolation  for  the  species  mass  fractions  is  not  required. 


SECTION  IV 


SOLUTION  OF  THE  FLUID  DYNAMIC  EQUATIONS 


1.  GENERAL 

In  Section  III,  the  fluid  dynamic  equations  were  defined  as  those 
equations  which  govern  the  fluid  dynamic  variables  u,  v,  P,  and  p. 

These  equations  are  the  global  continuity  equation,  the  component 
momentum  equations,  and  the  energy  equation.  Recall  from  Section  III 
that  in  the  inconsistent  scheme,  the  solutions  to  the  fluid  dynamic 
equations  are  advanced  one  time  step  while  the  distributions  of  the 
species  mass  fractions  throughout  the  nozzle  are  held  constant.  This 
section  describes  specifically  how  the  solutions  to  the  fluid  dynamic 
equations  are  determined. 

The  mesh  points  in  the  computational  plane  may  be  categorized  into 
four  groups:  interior,  inlet,  wall,  and  exit  points  (see  Figure  4). 

The  inlet,  wall,  and  exit  points  are  know*,  collectively  as  boundary 
points.  In  the  following  discussion,  th<»  numerical  treatment  used  for 
each  type  of  mesh  point  is  presented. 

2.  INTERIOR  POINTS 

Interior  mesh  points  are  computed  using  MacCormack's  method 
(Ref.  29);  an  explicit,  second-order  accurate,  two-step  finite  dif¬ 
ference  scheme.  The  fluid  dynamic  equations  are  employed  in  non¬ 
conservation  form.  Moretti  (30)  showed  that  using  the  conservation 
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Interior  (No 


computational  mesh  points 


form  decreased  computational  efficiency.  Backward  differences  are 
used  on  the  first  step  and  forward  differences  are  used  on  the  second 


step. 


As  an  example,  consider  the  partial  differential  equation 


u .  =  -uu 
t  x 


where  u  =  u(x,t) 


Application  of  MacCormack 1 s  method  yields  the  following  finite  differ 
ence  equations: 


Stepl:  *u5-&UX-u"-l> 


Step  2: 


[,/!  .  At  -n+l  (-n+1  .  flti+1  j  +0V+']/2 
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where  i  denotes  mesh  points  along  x,  n  denotes  the  time  step,  and  the 
tilde  denotes  values  calculated  by  the  first  step.  Similarly,  appli¬ 
cation  of  MacCormack's  method  to  equation  (7)  for  planar  flow  yields 
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where  L  and  M  denote  axial  and  radial  mesh  points,  respectively,  N 
denotes  the  time  step,  and  the  tilde  denotes  values  calculated  by  the 
first  step. 

When  the  nozzle  does  not  have  a  centerbody,  the  centerline  points 
are  treated  as  interior  points  and  solutions  at  these  points  are  com¬ 
puted  by  requiring  symmetry  about  the  centerline.  The  velocity  com¬ 
ponent  v  and  derivatives  of  u,  P,  and  p  with  respect  to  n  are  all  set 
to  zero.  The  derivative  of  v  with  respect  to  n  is  evaluated  by  a 
one-sided,  second-order  accurate  formula. 

3.  BOUNDARY  POINTS 

All  boundary  points  in  the  present  research  are  computed  using  a 
second-order  accurate  reference-plane  characteristic  scheme.  A 
reference-plane  scheme  is  used  rather  than  a  bicharacteristic  scheme 
because  the  increased  complexity  and  computational  time  of  the  bichar¬ 
acteristic  scheme  is  not  warranted  in  view  of  the  accuracy  limitations 
associated  with  the  chemical  kinetics  model.  In  reference-plane  char¬ 
acteristic  schemes,  derivatives  with  respect  to  one  of  the  independent 
variables  are  approximated  and  treated  as  forcing  terms,  thus  reducing 
the  number  of  independent  variables  in  the  problem  by  one.  For  example, 
in  the  constant  n  reference- plane  scheme,  all  derivatives  with  respect 
to  n  are  approximated  (by  MacCormack's  method)  and  are  placed  on  the 
right-hand  side  of  the  equal  sign  in  equations  (13)  to  (16).  Charac¬ 
teristic  relations  are  then  derived  for  the  resulting  equations.  A 
constant  n  reference-plane  method  is  used  for  inlet  and  exit  points 
while  a  constant  c  reference-plane  method  is  used  for  the  wall  points. 


The  characteristic  relations  for  these  reference-plane  methods  are 
derived  in  Appendix  C.  Figures  5  and  6  summarize  the  results  of  that 
Appendix. 

Cline  (4,  23)  computes  exit  points  by  linear  extrapolation  from 
the  first  and  second  adjacent  interior  points.  In  the  subject  problem, 
the  computational  mesh  extends  only  a  short  distance  into  the  super¬ 
sonic  flow  field  (just  far  enough  so  that  an  accurate  initial- 
data  surface  for  a  steady  flow,  method-of -characteristics  scheme  can 
be  computed).  The  constant  n  reference-plane  scheme  is  used  at  the 
exit  because  it  has  been  found  to  be  more  accurate  than  extrapolation 
in  the  "just  supersonic"  conditions  at  the  exit. 

In  the  following  paragraphs,  the  boundary  conditions  and  unit 
processes  for  inlet,  wall,  and  exit  points  are  described  in  turn. 

Inlet  Points 

Figure  5  illustrates  the  characteristic  relations  which  apply  for 
subsonic  flow  at  the  nozzle  inlet.  Note  that  in  this  case,  only  one 
characteristic  curve  is  contained  within  the  flow  field.  The  com¬ 
patibility  and  characteristic  equations  which  apply  to  the  subsonic 
inlet  case  are  (see  Appendix  C) 
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dP  -  pafdu  =  (<J>4  +  af  if/j  -  pa^Jdi 

(31) 

along  dc  =  (u-a^)dx 
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compatibility  equations 


characteristic  curves 


Figure  5.  Constant  n  reference- plane  characteristic  relations. 
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(34) 


♦4  *  -vP„  +  af2ypn  +  *k 

Since  there  are  4+n  unknowns  at  the  inlet  (u,  v,  P,  p,  and  C^, 
i=l,...,n)  and  only  one  compatibility  equation,  3+n  additional  condi¬ 
tions  must  be  specified. 

The  "correct"  specification  of  steady  flow  inlet  boundary  condi¬ 
tions  is  a  difficult  problem;  a  generally  preferred  choice  is  not 
known  at  the  present  time.  Of  the  several  different  sets  of  inlet 
boundary  conditions  which  have  been  tested  during  this  research, 
specification  of  0,  HQ,  P  ,  and  appears  to  give  the  best  results. 
Here,  e  is  the  inlet  flow  angle,  H0  is  the  total  enthalpy  (h  +  V2/2), 
PQ  is  the  total  pressure  based  on  frozen  composition  at  the  combustor 
exit,  and  are  the  species  mass  fractions.  The  specification  of 
total  conditions  at  the  inlet  places  an  upper  limit  on  the  energy  of 
the  flow  and  still  allows  flexibility  in  the  solution  for  the  static 
properties  at  the  inlet.  The  inlet  flow  angle  is  determined  by  exper¬ 
iment  or  from  the  following  procedure: 

1.  Determine  representative  values  of  the  ratio  of  frozen  speci¬ 
fic  heats  Yf  and  the  gas  constant  at  the  combustor  exit. 

2.  Generate  a  "long  inlet"  geometry  by  adding  six  to  ten  mesh 
points  upstream  of  the  nozzle  inlet,  simulating  a  constant- 
area  duct. 

3.  Using  the  program  of  the  subject  research  for  isentropic 
constant  specific  heat  ratio  flow,  compute  the  solution  to 
the  long  inlet  problem  described  in  the  preceding  two  steps. 
Inlet  boundary  conditions  for  the  long  inlet  problem  are: 


32 


a  uniform  distribution  of  zero  flow  angle,  and  total  pressure 
PQ  and  total  temperature  Tq  based  on  combustor  exit  conditions. 

4.  The  flow  angles  at  the  nozzle  inlet  for  the  converged  long 
inlet  problem  can  be  used  as  the  inlet  flow  angles  for  the 
reacting  flow  problem. 

Note  that  the  arbitrary  specification  of  zero  flow  angle  at  the  nozzle 
inlet  has  been  found  to  produce  unreasonable  velocity  profiles  at  the 
nozzle  inlet  for  some  nozzle  geometries. 

Examination  of  equation  (3)  shows  that  for  steady  flow,  the  total 

p 

enthalpy  (h  +  V  /2)  is  constant  along  streamlines.  Therefore,  the  dis¬ 
tribution  of  total  enthalpy  at  the  inlet  is  required  to  match  that  at 
the  combustor  exit. 

A  uniform  distribution  of  static  pressure  at  the  combustor  exit 
is  one  of  the  results  expected  from  the  combustor  analysis.  The 
specification  of  0,  HQ,  ,  and  uniform  static  pressure  was  one  of  the 
inlet  boundary  condition  sets  tested  during  this  research.  For  the 
cases  tested,  the  velocity  distributions  at  the  inlet  were  unreasonable 
with  very  low  (or  negative)  velocities  at  the  inlet  wall  point.  The 
use  of  total  pressure,  determined  from  conditions  at  the  combustor  exit 
and  based  on  stagnation  with  frozen  composition,  has  proved  effective 
in  computing  the  inlet  points. 

The  species  mass  fractions  are  specified  at  the  inlet.  Their 
values  are  determined  by  the  relaxation  of  the  species  continuity 
equations  between  combustor  exit  conditions  and  the  static  conditions 
at  the  inlet.  This  relaxation  process  is  described  in  detail  in  Sec¬ 
tion  V  and  Appendix  E. 
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For  the  case  of  supersonic  flow  at  the  inlet,  values  of  velocity, 
pressure,  density,  and  the  species  mass  fractions  are  all  specified 
since  downstream  conditions  do  not  propogate  upstream  in  a  supersonic 
flow. 

The  unit  process  at  the  inlet  for  subsonic  flow  is  based  on  a 
two-step,  predictor-corrector  method.  The  compatibility  equation  (31) 
and  the  specified  inlet  boundary  conditions  are  used  together  in  an 
iterative  procedure  to  compute  values  of  the  fluid  dynamic  variables 
at  the  solution  point.  The  iterative  procedure  is  described  first. 

Recall  that  0,  HQ,  P  ,  and  are  all  known  at  the  solution  point 
(point  3  in  Figure  5).  To  begin  the  iterative  procedure,  the  static 
temperature  T3  is  assumed.  Then  the  following  sequence  of  calculations 
is  performed  to  calculate  the  values  of  the  fluid  dynamic  variables: 

1.  Given  the  values  of  stavic  temperature  and  the  species  mass 
fractions  at  the  solution  point,  compute  the  gas  constant  R3, 
the  enthalpy  h3,  the  frozen  specific  heat  at  constant  pressure 
Cp3,  and  the  ratio  of  frozen  specific  heats  y^. 

2.  Using  the  definition  of  total  enthalpy  and  its  specified 
value  at  the  inlet,  compute  the  velocity  magnitude  at  the 
solution  point 

Q3  =  [2(Ho3  -  h3)]*  (35) 

Then  compute  the  frozen  speed  of  sound  a,  and  the  Mach  number 

r3 

at  the  solution  point. 

3.  Use  the  relation  between  total  and  static  pressure  at  the 
solution  point  to  compute  P3: 
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P3  -  Po3/[l  + 


(36) 


l-X^~  M32]Yf3-'<''f3-1> 

Now,  use  Pg,  Rg  and  Tg  in  the  thermal  equation  of  state  to 
compute  pg. 

4.  Compute  Ug  and  Vg  from  the  velocity  magnitude  at  the  solution 
point  Qg,  and  the  given  flow  angle. 

5.  Use  the  compatibility  equation  as  a  check  on  the  results  of 
the  preceeding  steps  [i.e.,  compute  using  the  compatibility 
equation  (31)  and  compare  with  Pg  computed  in  step  3  above]. 

6.  Iterate  on  Tg  using  the  secant  method  until  P^  computed  in 
steps  3  and  5  agree  to  within  a  specified  tolerance. 

This  iterative  procedure  is  used  as  part  of  each  predictor  and  corrector 
step  in  the  overall  unit  process  at  the  inlet.  The  predictor-corrector 
steps  are  described  next. 

In  order  to  use  the  compatibility  equation,  as  described  in  step 
5  of  the  iterative  procedure,  the  characteristic  curve  must  be  con¬ 
structed  and  the  compatibility  equation  must  be  written  in  finite 
difference  form.  This  is  done  by  replacing  the  differentials  in  equa¬ 
tion  (31)  with  differences  along  the  characteristic  curve.  In  the 
predictor  step,  all  coefficients  and  derivatives  are  evaluated  in  the 
initial-value  plane.  Given  the  time  step  and  the  estimate  for  (u-af), 
the  intersection  of  the  characteristic  curve  with  the  initial-value 
line  (in  the  constant  n  plane)  is  determined  using  the  finite  differ¬ 
ence  form  of  the  characteristic  curve  equation.  This  intersection  is 
point  2  in  Figure  5.  All  coefficients  in  equations  (31)  to  (34)  are 


evaluated  at  the  intersection  point  by  linear  interpolation  between 
points  A  and  B  in  Figure  5.  In  the  predictor  step,  the  derivatives 
with  respect  to  n  (in  the  ip  terms)  are  evaluated  using  backward  differ¬ 
ences  in  the  initial -value  plane.  When  all  coefficients  and  deriva¬ 
tives  are  evaluated  in  equations  (31)  to  (34),  the  iterative  procedure 
is  accomplished  and  the  predictor  step  is  complete.  This  is  done  for 
all  inlet  points  before  the  corrector  step  is  started. 

In  the  corrector  step,  the  characteristic  curve  must  be  con¬ 
structed  again.  Now,  the  coefficient  (u-af)  is  an  average  of  the 
values  at  the  solution  point  (from  the  predictor  step)  and  values  at 
point  2.  Also,  all  coefficients  in  equation  (31)  are  computed  as  the 
averages  of  values  at  the  solution  and  intersection  points.  Values 
at  point  2  are  again  determined  by  linear  interpolation.  Derivatives 
with  respect  to  n  in  the  ip  terms  are  evaluated  by  forward  differences 
at  the  solution  point  and  are  averaged  with  the  backward  difference 
approximations  in  the  initial-value  plane.  When  all  coefficients  and 
derivatives  in  equations  (31)  to  (34)  have  been  determined  by  the 
averaging  process  described  above,  the  iterative  procedure  is  accom¬ 
plished  and  the  overall  unit  process  at  the  inlet  is  complete. 

Wall  and  Centerbody  Points 

The  wall  and  centerbody  mesh  points  are  computed  using  the  con¬ 
stant  C  reference-plane  scheme  described  in  Appendix  C.  For  this 
scheme,  all  derivatives  with  respect  to  c  in  equations  (13)  to  (16) 
are  placed  on  the  right-hand  side  of  the  equal  sign  and  treated  as 
forcing  functions.  The  characteristic  and  compatibility  equations 


are  summarized  in  Figure  6. 

The  wall  slope  is  specified  at  each  mesh  point  along  the  wall 
(or  centerbody)  for  the  given  flow  geometry.  This  slope  provides  the 
single  boundary  condition  necessary  for  the  solution  of  the  fluid 
dynamic  equations  at  the  wall  (centerbody)  mesh  points.  For  example, 
consider  a  wall  point.  From  Figure  6,  the  appropriate  characteristic 
and  compatibility  equations  are 


6du  -  adv  =  (04*2  ~  “^Jd^ 
2 

dP  -  dp  =  ijj^dx 


along  dn  =  vdx 


dP 


paaf 

+  du 
a* 


p6af  o  P<*af 

+  dv  =  ^4  +  af  + 


p3a^. 


(37) 

(38) 

(39) 


v  =  u  tan  e 


along  dn  =  (v  +  afa*) 


(40) 


These  four  equations  in  the  four  fluid  dynamic  variables  are  written 
in  finite  difference  form  and  are  solved  in  a  predictor-corrector  pro¬ 
cedure  like  that  for  the  inlet  points.  Specifically,  equation  (40)  is 
substituted  into  equation  (37),  which  is  then  solved  for  u3  at  the 
solution  point.  Then  v3  is  obtained  from  equation  (40)  and  P3  is  com¬ 
puted  using  equation  (39).  Finally,  equation  (38)  is  solved  for  P3. 

Exit  Points 

Exit  points  are  computed  for  both  subsonic  and  supersonic  flow  by 
the  constant  n  reference-plane  procedure  (see  Figure  7).  For  subsonic 
flow,  the  exit  pressure  is  given  and  is  equal  to  the  ambient  pressure. 
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Subsonic 


C 

Supersonic 


Figure  7.  Subsonic  and  supersonic  exit  points 


The  characteristic  relations  in  this  case  are 

dv  =  if^dx  (41) 

along  dc  =  udx 

dP  -  af2dp  =  tp4dx  (42) 

dP  +  pafdu  =  (i|>4  +  +  Paf4'2)dx  (43) 

along  dc  =  (u  +  af)dx 

Equation  (41)  is  solved  for  Vg»  equation  (42)  is  solved  for  P3,  and 
equation  (43)  is  solved  for  u^. 

For  supersonic  flow  at  the  exit,  all  characteristic  curves  are 
within  the  flow  field  and  so,  another  compatibility  equation  is  added 
to  equations  (41)  to  (43)  above. 

o 

dP  -  pafdu  =  (<|>4  +  af  -  paf  ^2)dx 

(44) 

along  dc  =  (u  -  a^)dx 

Equations  (41)  to  (44)  are  four  equations  in  the  four  fluid  dynamic 
variables.  Equation  (41)  is  solved  for  v3  and  equations  (42)  to  (44) 
are  solved  iteratively  for  u^,  P^,  and  p3- 

For  both  subsonic  and  supersonic  flow,  a  predictor-corrector 
procedure,  analogous  to  that  for  the  inlet  mesh  points,  is  employed. 
Coefficients  and  derivatives  in  the  characteristic  relations  are  eval 
uated  in  the  initial -value  plane  for  the  predictor  step.  Derivatives 
in  the  \p  terms  are  approximated  by  backward  differences.  In  the 
corrector  step,  coefficients  and  derivatives  are  evaluated  in  the 
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SECTION  V 


SOLUTION  OF  THE  SPECIES  CONTINUITY  EQUATIONS 

1.  GENERAL 

Recall  from  Section  III  that  in  the  inconsistent  scheme,  the 
solutions  to  the  fluid  dynamic  equations  are  advanced  one  time  step 
while  the  distributions  of  the  species  mass  fractions  throughout  the 
nozzle  are  held  constant.  Section  IV  provides  the  details  of  this 
part  of  the  overall  numerical  algorithm.  Then,  the  new  flow  field  is 
assumed  to  be  steady  and  the  species  continuity  equations  are  inte¬ 
grated  along  as  many  as  four  streamlines  from  the  nozzle  inlet  to  the 
exit.  This  section  describes  specifically  how  the  integration  of  the 
species  continuity  equations  is  accomplished. 

It  is  well  known  that  integration  of  the  species  continuity  equa¬ 
tions  for  flows  near  chemical  equilibrium  requires  special  care  because 
the  equations  become  quite  "stiff"  (Refs.  3,  5,  and  31).  The  concept 
of  a  stiff  differential  equation  and  a  physical  explanation  for  the 
stiffness  of  the  species  continuity  equations  are  provided  in  Appendix 
D.  Many  numerical  techniques  have  been  proposed  for  the  solution  of 
the  species  continuity  equations  in  near  equilibrium  flows  (Refs.  32 
and  33).  Cline  and  Hoffman  (3)  analyzed  and  tested  a  number  of  pro¬ 
posed  schemes  in  their  analysis  of  three-dimensional,  steady,  non¬ 
equilibrium  f«ow.  They  concluded  that  explicit  schemes  and 


predictor-corrector  methods  are  stable  only  for  very  sma l  step  sizes. 
Also,  only  some  of  the  implicit  methods  they  tested  gave  adequate 
results.  A  method  proposed  by  Lomax  and  Bailey  (34)  was  chosen  as  the 
preferred  method  and  it  is  this  scheme  which  has  been  used  in  the 
subject  research.  The  method  is  second-order  accurate,  implicit,  and 
is  based  on  a  Taylor  series  expansion  of  the  species  continuity  equa¬ 
tions. 


2.  SECOND-ORDER,  IMPLICIT,  TAYLOR  EXPANSION 

The  second-order,  implicit  Taylor  expansion  method,  as  applied 
to  the  species  continuity  equations,  is  derived  in  Appendix  D.  The 
resulting  finite  difference  equation  is 


cn+i =  cu  +  2v£  [fu(3'  vA1)+  ap*  (p*+rp*)  +  aT*  (WV 


n  3f . 


(45) 
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(C  -  C  ■  )]  (i=l,...,n) 


jii  3Cj  '  «+1  J* 


where  h  =  s^+1  -  s^,  s  is  the  position  along  the  streamline,  fi  =  a. j/p, 
V  is  the  velocity  magnitude,  and  the  subscripts  i  and  £+1  denote 
points  along  the  streamline.  The  partial  derivatives  in  equation  (45) 
are  determined  analytically.  When  expanded,  equation  (45)  becomes  a 
system  of  n  simultaneous,  linear,  algebraic  equations  in  the  n  unknowns 
=  C^  +j-C^.  This  system  of  equations  is  solved  using  Gauss  eli¬ 
mination  with  scaling  and  partial  pivoting. 

For  flows  which  are  not  near  chemical  equilibrium,  a  modified 
Euler  method  is  available  for  integration  of  the  species  continuity 
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equations.  The  adequacy  of  this  explicit  integration  method  must  be 
determined  on  a  trial-and-error  basis. 

3.  STREAMLINE  TRACING 

The  species  continuity  equations  are  integrated  along  streamlines 
in  the  inconsistent  scheme.  Since  streamlines  will  not  pass  through 
computational  mesh  points  in  general  (except  along  the  wall  and  cen- 
terbody/centerl ine) ,  it  is  necessary  to  trace  streamlines  through  the 
flow  field  before  the  species  continuity  equations  are  integrated.  The 
computer  program  developed  during  the  subject  research  allows  the  user 
to  trace  as  many  as  four  streamlines  through  the  flow  field.  The  two 
streamlines  along  the  wall  and  centerbody/centerl ine  must  be  used  and 
one  or  two  more  can  be  used  with  origins  at  arbitrary  mesh  points  along 
the  nozzle  inlet.  The  use  of  more  than  four  streamlines  in  the  over¬ 
all  algorithm  would  lead  to  prohibitively  long  computational  times. 

Figure  8  illustrates  the  streamline  tracing  procedure.  A 
predictor-corrector  method  is  used.  Starting  at  the  point  (1,M,N) 
at  the  nozzle  inlet,  the  values  of  u  and  v  at  this  point  are  used  to 
project  a  streamline  to  the  intersection  with  the  line  of  mesh  points 
at  L=2.  Then,  values  of  (u*^*  v*2^p  at  ^e  intersection  point  are 
found  by  linear  interpolation  between  known  values  at  the  adjacent 
mesh  points  at  L=2.  In  the  corrector  step,  the  values  of  (u*£»  v*2^p 
are  averaged  with  u  and  v  at  (1,M,N)  and  the  streamline  is  again  pro¬ 
jected  to  the  intersection  at  L=2.  Linear  interpolation  is  used  to 
find  (u*2>  v*2)c*  This  predictor-corrector  method  is  applied  at  each 
L  station  until  the  streamline  has  been  traced  to  the  nozzle  exit. 


4.  ERROR  CONTROL 

The  implicit  scheme  described  previously  has  been  found  to  be 
stable  for  even  large  step  sizes  in  both  Cline  and  Hoffman's  work  (3) 
and  in  the  subject  research.  However,  the  method  is  subject  to 
truncation  error,  and  thus  consideration  must  be  given  to  step  size 
control . 

In  the  derivation  of  equation  (45),  third-  and  higher-order  terms 
have  been  neglected.  The  objective  of  the  error  control  scheme  is  to 
estimate  the  third-order  term  and  (by  adjusting  the  step  size)  to 
insure  that  the  ratio  of  this  term  to  the  computed  term  [the  right- 
hand  side  of  equation  (45)]  is  less  than  a  specified  tolerance.  An 
expression  for  this  ratio  is  derived  in  Appendix  D  and  is  restated 
here: 

K,nJ,  -  2K.„  +  K.„  , 

RATIO.  =  (i=l . n)  (46) 

1  6Kifc 

Note  that  RATIO  has  distinct  values  for  each  of  the  chemical  species 
in  the  gas  mixture  at  each  step  in  the  integration.  Also,  RATIO  can 
be  computed  only  after  three  integration  steps  have  been  taken. 

The  error  control  scheme  is  implemented  by  placing  intermediate 
points  between  grid  points  along  the  streamlines  (see  Figure  9).  An 
estimate  is  made  of  the  number  of  intermediate  points  (NINT)  required 
between  L*  and  L+l*.  Values  of  temperature,  density,  and  velocity  at 
L*  and  L+l*  are  determined  by  linear  interpolation  between  known 
values  at  adjacent  mesh  points  at  each  L  station.  The  values  of  these 
same  variables  at  the  intermediate  points  are  determined  by  linear 
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interpolation  between  the  corresponding  values  at  L*  and  L+l*.  Then, 
the  species  continuity  equations  are  integrated  through  the  inter¬ 
mediate  points  from  L*  to  L+l*  and  RATIO  is  computed  for  each  of  the 
chemical  species.  If,  after  the  last  intermediate  integration  step, 
the  maximum  value  of  RATIO  is  less  than  a  user  specified  tolerance, 
then  the  integration  proceeds.  If  the  required  tolerance  is  not 
achieved,  then  the  number  of  intermediate  points  is  doubled  and  the 
integration  procedure  is  restarted  from  L*.  Also,  if  the  maximum 
value  of  RATIO  is  two  orders  of  magnitude  less  than  the  specified 
tolerance,  then  the  number  of  intermediate  points  is  halved  before  the 
integration  from  L+l*  to  L+2*  proceeds. 

Numerical  experiments  have  been  performed  to  illustrate  the  use  of 
the  error  control  scheme.  Oetails  are  provided  in  Appendix  D  and  the 
results  are  summarized  here. 

1.  Most  of  the  computational  time  required  to  solve  the  subject 
problem  is  related  to  integration  of  the  species  continuity 
equations.  Therefore,  the  total  execution  time  is  nearly 
directly  proportional  to  the  number  of  intermediate  points 
(NINT). 

2.  Failure  to  use  intermediate  points  will,  in  general,  yield 
poor  results  for  species  mass  fractions  and  the  energy 
source  term  4^. 

3.  Integrating  from  equilibrium  conditions  with  very  small 
gradients  in  temperature  and  density  and  a  small  value  of  the 
specified  tolerance  on  RATIO  can  yield  very  large  numbers  of 
intermediate  points. 


4.  For  equilibrium  initial  conditions,  NINT  is  not  particularly 
sensitive  to  the  sign  or  magnitude  of  temperature,  density, 
or  velocity  gradients  between  L*  and  L+l*  as  long  as  the 
conditions  in  conclusion  (3)  do  not  occur. 

5.  For  nonequilibrium  initial  conditions,  NINT  is  sensitive  to 
gradients  in  temperature  and  velocity  while  it  is  only  mildly 
sensitive  to  the  density  gradient  between  L*  and  L+l*.  If 
the  property  gradients  between  these  points  do  not  match  the 
corresponding  gradients  between  L-l*  and  L*,  more  intermediate 
points  are  likely  to  be  required  to  achieve  a  specified 
tolerance. 

6.  Conclusions  1  to  5  appear  to  be  valid  for  the  several  chemis¬ 
try  systems  investigated  in  this  research. 

Because  of  conclusion  3,  an  upper  limit  of  20  intermediate  points 
is  used  in  the  computer  program  for  the  subject  problem.  Experience 
indicates  that  a  relatively  large  tolerance  (5  to  10%)  should  be  used 
until  the  solution  is  near  convergence  and  then  the  tolerance  can  be 
reduced  to  the  desired  level.  Also,  the  user  has  the  option  of  speci¬ 
fying  the  number  of  intermediate  points  thereby  overriding  the  error 
control  scheme. 

5.  RELAXATION  AT  THE  INLET 

The  choice  of  inlet  boundary  conditions  was  described  in  Section 
IV.  Recall  that  the  flow  angle  e,  the  total  enthalpy  HQ,  the  total 
pressure  based  on  combustor  exit  conditions  and  frozen  composition  P  , 
and  the  species  mass  fractions  are  all  specified  at  each  inlet  mesh 


point.  The  use  of  total  conditions  as  inlet  boundary  conditions 
establishes  the  energy  level  of  the  flow  and  allows  some  flexibility 
in  the  distribution  of  static  properties  across  the  inlet.  As  a 
result,  the  distributions  of  the  static  properties  across  the  inlet 
will  not  necessarily  match  those  across  the  combustor  exit.  Discon¬ 
tinuities  in  static  properties  between  the  combustor  exit  and  the 
nozzle  inlet  create  the  need  for  a  relaxation  of  the  species  mass 
fractions  between  these  two  points.  A  detailed  description  of  the 
relaxation  process  is  provided  in  Appendix  E.  Here,  the  justifica¬ 
tion  for  the  relaxation  and  the  relaxation  process  are  summarized. 

Consider  a  particle  (a  small  mass  of  the  reactive  mixture  of 
gases)  at  the  combustor  exit.  Also,  note  that  in  functional  form 
(see  Appendix  A) 

oi  =  a.(p,T,Ci)  (i=l . n)  (47) 

=  ^k^T,Gi’Ci^  =  lMp,T’{V  (48) 

Since  all  the  static  properties  and  the  species  mass  fractions  are 
known  at  the  location  of  the  particle,  the  species  source  functions 
a..  and  the  energy  source  term  ^  can  be  computed  at  this  point.  In 
a  real  flow,  there  can  be  no  discontinuities  in  the  streamwise  dis¬ 
tributions  of  the  static  properties  and  species  mass  fractions  at  the 
junction  between  the  combustor  exit  and  the  nozzle  inlet.  Therefore, 
from  equations  (47)  and  (48),  the  distributions  of  o..  and  along  the 
flow  direction  must  also  be  continuous  at  this  junction.  The  dis¬ 
continuities  in  the  static  properties  described  in  the  preceding 
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paragraph  are  a  consequence  of  the  problem  model;  that  is,  an  attempt 
to  link  a  two-dimensional  nozzle  analysis  to  a  quasi-one-dimensional 
combustor  analysis.  Since  the  discontinuities  do  not  represent  reality, 
their  effect  should  be  manifest  at  the  inlet  before  the  species  con¬ 
tinuity  equations  are  integrated  through  the  flow  field.  This  is 
accomplished  by  allowing  the  species  mass  fractions  at  the  combustor 
exit  to  adjust  to  the  static  conditions  at  the  nozzle  inlet.  The 
species  continuity  equations  are  relaxed  between  conditions  at  the 
combustor  exit  and  those  at  the  inlet  for  a  period  of  time  that  mini¬ 
mizes  the  change  in  between  these  two  points.  In  effect,  a  short 
distance  is  spliced  into  the  junction  between  the  combustor  exit  and 
the  nozzle  inlet  so  that  the  discontinuities  in  static  properties 
produced  by  the  mathematical  model  can  be  replaced  with  linear  grad¬ 
ients.  The  distance  is  chosen  so  that  the  chemical  nature  of  the  flow 
at  the  combustor  exit  (as  represented  by  the  energy  source  term  i|^) 
is  preserved  at  the  nozzle  inlet  to  the  extent  possible.  Figure  10 
illustrates  the  effect  of  relaxation  time  on  ^  for  an  H-F  propellant 
system.  The  initial  conditions  are  near  equilibrium  and  curves  for 
several  different  temperature  gradients  are  shown.  Figure  11  shows 
relaxation  curves  for  the  same  propellant  system  and  gradients  but 
with  nonequilibrium  initial  data.  Several  important  conclusions 
regarding  relaxation  curves  from  Appendix  E  are  restated  here: 

1.  All  curves  of  versus  relaxation  time  approach  zero  as 
relaxation  time  increases.  This  is  because  the  chemical 
systems  move  toward  an  equilibrium  condition  (where  the 
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species  source  functions  go  to  zero)  with  increasing  relaxa¬ 
tion  time. 

2.  At  high  temperatures  and  near  equilibrium  conditions,  ^  is 
quite  sensitive  to  variations  in  temperature  and  density. 

For  example,  a  20°R  temperature  perturbation,  with  no  relaxa¬ 
tion,  can  change  the  computed  value  of  ^  by  several  orders 
of  magnitude  (see  Figure  10). 

3.  For  the  chemistry  systems  and  static  property  gradients 
studied  in  this  research,  relatively  small  changes  in  composi¬ 
tion  occur  during  the  relaxation  from  combustor  exit  to 
nozzle  inlet  conditions.  Typically,  species  mass  fractions 
change  by  less  than  two  percent. 

Failure  to  relax  the  species  mass  fractions  at  the  inlet  leads  to 
numerical  difficulties  for  two  primary  reasons.  First,  without  relax¬ 
ation,  the  data  set  (p.T.C^ ,  i=l,...,n)  at  the  inlet  is  inconsistent. 

As  noted  above,  the  energy  source  term  ^  (and  species  source  function) 
can  then  differ  by  several  orders  of  magnitude  from  its  value  at  the 
combustor  exit.  This  causes  a  "hard  start"  condition  in  integration  of 
the  species  continuity  equations  away  from  the  inlet.  Second,  if 
the  species  continuity  equations  are  not  relaxed  before  the  inlet, 
they  will  relax  within  the  computed  flow  field.  This  can  produce  dis¬ 
tributions  of  within  the  flow  with  dramatic  oscillations  near  the 
inlet;  this  ij^  distribution,  in  turn,  distorts  the  flow  field. 

The  proper  relaxation  time  is  determined  as  part  of  the  overall 
numerical  algorithm.  The  scheme  is  illustrated  in  Figure  12  which 


Figure  12.  Calculate 
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shows  a  typical  4^  versus  relaxation  time  curve.  In  that  figure,  4^ 
is  computed  for  the  conditions  at  the  combustor  exit,  4^  is  computed 
for  the  species  mass  fractions  at  the  combustor  exit  and  density  and 
temperature  at  the  nozzle  inlet  with  zero  relaxation  time,  and  4^ 
is  computed  for  the  same  conditions  as  4^  but  with  a  relaxation  time 
At.  The  values  of  4^  corresponding  to  each  inlet  mesh  point  are 
computed  as  part  of  the  initialization  procedure  for  the  overall  algo¬ 
rithm.  They  are  only  computed  once  and  then  are  stored.  Consider 
time  level  N  in  the  solution  of  the  subject  problem.  4kQ  is  computed 
at  each  streamline  origin  for  the  static  conditions  prevailing  at  the 
inlet.  There  is  a  minimum  relaxation  time  and  if  |4.qI  <  | > 
then  the  minimum  relaxation  time  is  used.  This  is  because  4^  is  known 
to  approach  zero  with  increasing  relaxation  time.  Otherwise,  the 
species  continuity  equations  are  relaxed  through  time  At1  (see  Figure 
12)  which  is  the  relaxation  time  computed  for  a  given  streamline 
during  previous  time  steps.  Then  the  species  continuity  equations  are 
integrated  along  the  given  streamline  to  the  nozzle  exit.  Ne^t,  a 
check  is  made: 


l^kAt  '  ^kR  I 
I*k0  ‘  ♦kR  1 


? 

<  .05 


(49) 


If  the  check  is  satisfied,  then  At^  is  not  changed  and  it  is  saved  for 
use  during  the  next  time  step  in  the  overall  algorithm.  However,  if 
the  check  is  not  satisfied,  then  a  straight  lint  :s  extended  from  4 
through  4^At  to  intersect  4'^.  The  relaxation  time  At^  corresponding 
to  that  intersection  is  used  during  the  next  time  step.  The  same 
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procedure  continues  for  subsequent  time  steps  as  shown  in  Figure  12 
until  the  check  is  satisfied. 


SECTION  VI 


OVERALL  NUMERICAL  ALGORITHM 


1.  GENERAL 

The  overall  numerical  algorithm  consists  of  the  repetitive  and 
alternate  application  of  the  procedures  described  in  Sections  IV  and 
V;  the  solutions  to  the  fluid  dynamic  equations  are  advanced  one  time 
step  for  a  given  distribution  of  the  species  mass  fractions  throughout 
the  nozzle,  and  then  the  flow  field  is  assumed  to  be  steady  while  the 
species  continuity  equations  are  integrated  along  as  many  as  four 
streamlines  from  the  nozzle  inlet  to  the  exit.  Note  that  for  some 
problems,  the  fluid  dynamic  equations  may  be  advanced  more  than  one 
time  step  before  the  species  continuity  equations  are  integrated.  If 
this  is  possible,  the  total  computational  time  may  be  reduced  signifi¬ 
cantly  (see  Section  VII).  As  discussed  in  Section  III,  this  overall 
scheme  is  inconsistent  in  time  in  the  treatment  of  the  species  con¬ 
tinuity  equations,  but  it  becomes  consistent  at  the  steady  state  limit. 
Also,  it  is  necessary  to  relax  the  species  mass  fractions  between  the 
specified  static  conditions  at  the  combustor  exit  and  those  at  the 
nozzle  inlet  before  integrating  the  species  continuity  equations 
through  the  flow.  Discontinuities  in  static  properties  at  the  junction 
between  the  combustor  exit  and  the  nozzle  inlet  occur  because  of  the 
choice  of  inlet  boundary  conditions  for  the  problem  model.  An  error 


control  scheme,  which  limits  the  size  of  the  truncation  error  in  the 
integration  of  the  species  continuity  equations,  is  available.  Con¬ 
vergence  is  achieved  when  the  computed  values  of  the  fluid  dynamic 
variables  and  the  species  mass  fractions  no  longer  change  with  time. 

Two  sets  of  mesh  points  are  used  in  the  overall  algorithm,  as 
illustrated  in  Figure  13.  During  the  first  part  of  the  algorithm,  the 
solutions  for  u,  v,  P,  and  p  are  advanced  one  time  step  at  all  compu¬ 
tational  mesh  points  in  a  predictor-corrector  procedure  (see  Section 
IV).  Then,  the  flow  is  assumed  to  be  steady,  and  streamlines  are 
traced  through  the  flow  field  to  locate  the  streamline  mesh  points. 

The  values  of  static  properties  at  the  streamline  mesh  points  which 
are  required  for  integration  of  the  species  continuity  equations  are 
determined  by  interpolation  between  adjacent  computational  mesh  points 
at  each  L  station.  In  addition  to  computing  the  distributions  of  the 
species  mass  fractions  along  streamlines  when  the  species  continuity 
equations  are  integrated,  the  energy  equation  source  term  4^.  the  ratio 
of  frozen  specific  heats  y^,  the  gas  constant  R,  and  the  total 
enthalpy  Hq  are  also  computed  and  stored  for  each  streamline  mesh 
point.  The  variables  y^,  R,  and  4^  are  all  needed  for  the  solution  of 
the  fluid  dynamic  equations.  The  total  enthalpy  distribution  is  used 
as  a  check  on  the  validity  and  accuracy  of  the  converged  solution. 

The  values  of  these  variables  at  the  computational  mesh  points  are 
determined  by  linear  interpolation  between  adjacent  streamline  mesh 
points  at  each  L  station. 
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mesh  points. 


2.  INITIAL -VALUE  SURFACE 

The  initial -value  surface  must  be  a  smooth  distribution  of  the 
fluid  dynamic  variables  u,  v,  P,  and  p  throughout  the  computational 
mesh.  It  can  be  read  in  or  computed  internally  within  the  computer 
program  for  the  subject  research.  Three  different  methods  have  been 
investigated  for  internal  generation  of  the  initial-value  surface: 

1.  Obtain  a  one-dimensional  solution  to  the  problem  of  interest  using 
the  program  of  Reference  (10).  Then,  use  the  results  of  this 
analysis  to  provide  the  axial  distributions  of  pressure,  density, 
temperature  and  velocity  magnitude  for  the  initial-value  surface. 
The  local  flow  angle  0  and  velocity  components  at  any  point  in  the 
initial-value  surface  are  estimated  from  the  given  geometry  (i.e., 
interpolate  between  wall  and  centerbody  values  for  0  at  each  L 
station) . 

2.  Obtain  a  one-dimensional,  isentropic,  constant  specific  heats 
solution  for  pressure,  density,  and  velocity  magnitude.  Values  of 
R  and  y  for  this  analysis  are  representative  values  at  the  combus¬ 
tor  exit.  The  local  flow  angle  and  velocity  components  are  esti¬ 
mated  as  in  method  1. 

3.  Using  R  and  Y  as  in  method  2,  compute  the  two-dimensional,  isen¬ 
tropic,  constant  specific  heats  solution  for  the  subject  problem 
and  use  it  as  the  initial -value  surface. 

For  all  three  methods,  the  distributions  of  species  mass  fractions  are 
determined  by  integrating  the  species  continuity  equations  through  the 
given  flow  field  from  the  known  conditions  at  the  nozzle  inlet. 
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Method  1  requires  the  use  of  two  computer  programs  and  the  trans¬ 
fer  of  data  from  one  to  the  other.  Also,  there  is  no  provision  in  the 
program  of  Reference  (10)  for  nonuniform  data  at  the  nozzle  inlet. 
Methods  2  and  3  can  be  accomplished  using  only  the  computer  program  of 
the  subject  research  and  can  be  used  for  any  set  of  data  at  the  nozzle 
inlet.  Method  2  is  simpler  to  implement  than  method  3  and,  since  there 
appears  to  be  little  difference  in  overall  computational  time  to  con¬ 
vergence  for  the  two  methods,  method  2  is  preferred.  All  three  methods 
have  yielded  essentially  the  same  converged  solutions  for  the  cases 
tested. 

3.  BOUNDARY  CONDITIONS 

The  choice  of  steady  flow  boundary  conditions  has  been  discussed 
previously  in  Section  IV.  Recall  that  for  subsonic  flow,  the  flow 
angle  6,  the  total  enthalpy  Hq,  the  total  pressure  based  on  combustor 
exit  conditions  PQ,  and  the  species  mass  fractions  are  specified 
at  the  inlet;  the  flow  angle  0  is  specified  along  the  wall  and  center- 
body;  and,  the  exit  pressure  is  set  equal  to  the  ambient  pressure. 

For  supersonic  flow,  all  static  conditions  and  species  mass  fractions 
are  specified  at  the  inlet  and  the  flow  angle  is  specified  along  the 
wall  and  centerbody. 

It  is  essential  that  the  data  at  the  combustor  exit  be  consistent 
with  the  chemical  kinetics  model  that  is  employed.  This  fact  has  been 
noted  by  Cline  and  Hoffman  (3)  and  has  been  established  again  in  the 
subject  research.  An  inconsistent  set  of  data  will  yield  computed 
values  of  the  species  source  functions  which  can  be  several  orders  of 


61 


magnitude  in  error.  As  a  result,  an  attempt  to  integrate  the  species 
continuity  equations  away  from  the  combustor  exit  yields  immediate 
numerical  difficulty. 

4.  STEP  SIZE  AND  STABILITY 

No  attempt  has  been  made  to  perform  a  stability  analysis  for  the 
overall  numerical  algorithm  as  applied  to  the  governing  equations  for 
two-dimensional,  unsteady,  nonequilibrium,  chemically  reacting  flow. 

The  stability  of  the  treatment  of  the  equations  governing  the  fluid 
dynamic  variables  and  the  stability  of  the  species  continuity  equation 
scheme  are  treated  separately.  The  stability  of  the  overall  numerical 
algorithm  has  been  verified  by  numerical  experiment. 

The  scheme  used  to  advance  the  solutions  for  the  fluid  dynamic 
variables  is  subject  to  the  CFL  stability  restriction.  This  requires 
that  the  finite  difference  domain  of  influence  must  be  at  least  as 
large  as  the  continuum  domain  of  influence.  It  ensures  that  the  speed 
of  propogation  of  numerical  disturbances  (truncation  error  for  example) 
everywhere  exceeds  the  speed  of  propogation  of  disturbances  in  the 
flow  (i.e.,  the  speed  of  sound  in  a  compressible  flow).  Application 
of  the  CFL  criterion  to  two-dimensional,  unsteady  flow  yields 


At  £ 


[<V  +  a)(-i+-ij)J] 
hC  Arf 


In  practice,  the  step  size  is  computed  as 


Ax  = 


[(V  +  aH-V-^)*] 

A<^  Arf 


(50) 


(51) 
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where  A  ranges  from  0.4  to  1.6  depending  on  the  geometry  of  the  flow 
[Reference  (23)].  Any  study  of  the  time  step  multiplier  A  should  be 
performed  for  a  nonreacting  flow  similar  to  the  flow  with  finite-rate 
kinetics.  The  proper  choice  of  A  can  accelerate  convergence.  The 
author's  experience  indicates  that  the  value  of  A  used  for  a  represen¬ 
tative  nonreacting  flow  will  also  work  for  the  reacting  flow  problem. 

The  second-order,  implicit  scheme  used  for  integration  of  the 
species  continuity  equations  has  been  found  to  be  stable  for  all  step 
sizes.  Note  that  the  step  size  for  integration  of  these  equations  is 
not  related  to  Ax  in  equation  (51)  since  the  species  continuity  equa¬ 
tions  are  integrated  along  streamlines  in  space  while  the  fluid  dynami 
equations  are  integrated  from  one  time  plane  to  the  next. 

5.  THERMOCHEMICAL  MODELS 

In  addition  to  the  chemically  reacting  flow  of  a  mixture  of  ther¬ 
mally  perfect  gases,  two  other  thermochemical  models  can  also  be 
analyzed  using  the  computer  program  of  the  subject  research.  These 
models  are: 

1.  Isentropic  flow  of  a  thermally  and  calorically  perfect  gas 
(constant  y). 

2.  Isentropic  flow  of  a  gas  whose  equation  of  state  is  input  in  tabu¬ 
lar  form.  Two-dimensional  frozen  and  equilibrium  solutions  can 

be  obtained  using  this  model. 
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SECTION  VII 


VERIFICATION  AND  RESULTS 


1.  VERIFICATION 

At  the  present  time,  there  are  no  other  methods  available  for  the 
analysis  of  two-dimensional,  nonequilibrium,  chemically  reacting  sub¬ 
sonic  and  transonic  nozzle  flows.  Therefore,  it  is  not  possible  to 
make  comparisons  with  existing  results.  However,  a  number  of  tests 
have  been  performed  to  verify  the  results  of  the  subject  research 
as  discussed  in  the  following  paragraphs. 

The  scheme  for  integration  of  the  species  continuity  equations 
must  generate  accurate  concentration  profiles  for  given  distributions 
of  the  fluid  dynamic  variables  along  the  nozzle  axis.  In  order  to 
test  this  capability  in  the  subject  computer  program,  one-dimensional 
analyses  of  several  kinetics  problems  were  performed  using  the  program 
of  Reference  (10).  These  problems  are  discussed  in  Appendix  E.  Then, 
using  the  distributions  of  p,  T,  and  V  from  the  results,  the  species 
continuity  equations  were  integrated  through  the  given  flow  fields  by 
the  technique  of  the  subject  research.  The  resulting  concentration 
profiles  were  compared  with  those  computed  by  Reference  (10).  For  the 
C-H-O-N  system,  an  "error"  of  0.35  percent  in  the  mass  fraction  of 
monatomic  hydrogen  at  the  nozzle  throat  constituted  the  most  signifi¬ 
cant  deviation  from  the  Reference  (10)  solution.  Applying  the  same 


procedure  to  the  H-F  system  yielded  “errors"  at  the  throat  which  were 
typically  less  than  two  percent.  The  worst  case  was  a  3.4  percent 
deviation  in  the  mass  fraction  of  monatomic  flourine.  Figure  14 
illustrates  the  species  mass  fraction  profiles  for  the  H-F  system  as 
computed  by  the  integration  technique  of  the  subject  research. 

There  are  four  results  which  must  be  true  of  a  converged  solution 
to  the  subject  problem  and  which  are  therefore  useful  in  determining 
the  validity  of  the  computed  results.  First,  the  law  of  conservation 
of  mass  requires  that  the  mass  flow  at  each  axial  position  along  the 
nozzle  must  be  constant.  Values  of  mass  flow  at  the  nozzle  inlet, 
throat,  and  exit  are  computed  and  made  available  as  part  of  the  program 
output.  Second,  in  the  steady  state  limit,  the  total  enthalpy  Hq 
must  be  constant  along  streamlines  in  the  flow.  A  "total  enthalpy 
error"  has  been  defined  so  that  this  fact  can  be  used  as  a  measure  of 
the  accuracy  of  the  solution. 

H  Error  i  (52) 

0  V  c  _  V  c 
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where  HQL  is  the  total  enthalpy  at  axial  position  L  along  a  given 
streamline,  Hq^  is  the  total  enthalpy  at  the  combustor  exit  for  the 
same  streamline,  and  VE  and  Vj  are  the  velocity  magnitudes  at  the 
nozzle  exit  and  inlet,  respectively,  along  the  streamline.  Profiles 
of  total  enthalpy  error  for  the  cases  studied  are  presented  with  the 
discussion  of  results.  Third,  the  sum  of  the  species  mass  fractions 
must  be  one  throughout  the  flow.  This  fact  is  not  used  explicitly 


I 


5 


10  15 

Axial  Position 


20 


Figure  14.  Concentration  profiles  for  the  H 
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F  system. 


in  the  problem  formulation  or  in  the  numerical  algorithm.  However,  the 
sum  of  the  species  mass  fractions  at  each  axial  mesh  point  along  the 
kinetics  streamlines  is  computed  and  output.  The  fourth  result  which 
must  be  true  of  the  converged  solution  is  that  the  finite-rate  kinetics 
results  must  lie  between  the  equilibrium  and  frozen  solutions  since 
those  cases  correspond  to  infinite  and  zero  reaction  rates,  respec¬ 
tively. 

As  a  further  test,  the  results  of  the  subject  technique  have  been 
mass-averaged  (in  the  radial  direction)  and  compared  with  an  accepted 
one-dimensional  solution  [Reference  (10)].  Mass-averaged  temperature 
profiles  are  compared  with  one-dimensional  temperature  profiles  in 
the  discussion  of  results. 

In  the  final  verification  step,  the  subject  technique  has  been 
applied  to  a  supersonic  problem  and  the  results  have  been  compared 
with  those  from  a  two-dimensional  method-of-characteri sties  solution 
[Reference  (2)].  The  results  of  an  H-F  system  subsonic- transonic 
analysis  were  used  to  generate  the  initial  data  line  for  both  the 
method-of -characteristics  solution  and  the  solution  by  the  subject 
method.  Eight  equally  spaced  mesh  points  were  used  along  the  initial 
data  line.  A  comparison  of  the  results  is  presented  in  Figure  15  (TDK 
corresponds  to  the  char' 'teristics  solution).  Note  the  excellent 
agreement  in  the  temperature  and  concentration  profiles  along  the 
wall  and  the  fair  agreement  of  the  centerline  temperature  profile. 
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Figure  15 
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.  Verification  with  the  TDK  computer  program 


2.  RESULTS 


Several  different  nozzle  geometries  and  chemistry  systems  have 
been  analyzed  during  this  research.  The  details  of  the  nozzle  geo¬ 
metries  and  reaction  mechanisms  for  the  systems  discussed  here  are 
presented  in  Appendix  E. 

C-H-O-N  System 

The  C-H-O-N  system  corresponds  to  a  realistic  application  of  the 
subject  technique  to  a  ramjet  nozzle  problem.  The  composition  of  the 
reactive  mixture  of  gases  at  the  nozzle  inlet  results  from  an  equili¬ 
brium  calculation  for  the  combustion  of  a  hydrocarbon  fuel  in  air. 

The  nozzle  size  and  shape  is  appropriate  for  a  ramjet  application.  It 
is  a  circular-arc  conical  nozzle  with  a  7.5  inch  inlet  radius,  «  4.75 
inch  throat  radius,  a  45  degree  angle  in  the  convergence,  and  a  15 
degree  angle  in  the  divergence.  There  are  12  chemical  specie,  and  8 
reactions  in  the  chemistry  model.  Figure  16  presents  the  results  of 
the  analysis  for  this  problem.  Note  that  the  mass-averaged  tempera¬ 
ture  for  the  finite-rate  kinetics  solution  does  fall  between  the 
equilibrium  and  f ozen  solutions  and  that  it  matches  the  ODK  solution 
[Reference  (10)]  quite  well.  Also  note  that  in  this  case,  there  is 
only  a  small  departure  from  the  frozen  solution.  For  both  the  isen- 
tropic  (y  =  constant)  and  kinetics  solutions,  a  mass  flow  difference 
of  0.8  percent  between  tne  inlet  and  the  throat  values  occurred  in 
the  converged  results.  For  the  kinetics  solution,  the  sum  of  the 
species  mass  fractions  deviaced  by  less  than  lx  10"^  along  both  the 
center  and  wall  streamlines.  The  profiles  of  total  enthalpy  error 


Temperature  profiles  for  the  C-H-O-N  system. 


are  presented  in  Figure  17.  Note  that  the  error  was  less  than  one 
percent  at  all  mesh  points  except  for  the  isentropic  solution  at  the 
last  wall  point.  This  particular  point  is  just  downstream  from  the 
junction  between  the  circular  arc  in  the  throat  region  and  the  conical 
section  in  the  nozzle  divergence.  The  author's  experience  indicates 
that  the  discontinuity  in  the  second  derivative  of  the  wall  coordi¬ 
nates  at  this  junction  can  cause  distortions  in  the  total  enthalpy 
profile.  Adding  three  more  mesh  points  in  the  axial  direction 
(L=16  to  L=18)  and  recomputing  the  isentropic  solution  reduces  the 
total  enthalpy  error  at  L=15  to  one  percent.  Note  in  Figure  17  that 
the  addition  of  finite-rate  chemical  kinetics  to  this  problem  has  not 
significantly  affected  the  Hq  profiles  with  respect  to  the  isentropic 
(y  =  constant)  results. 

With  18  mesh  points  along  the  axis,  two  chemical  kinetics  stream¬ 
lines,  and  five  intermediate  points  specified,  each  integration  of  the 
species  continuity  equations  for  the  C-H-O-N  system  required  approxi¬ 
mately  11  seconds  on  the  CDC  6500  computer.  The  addition  of  a  third 
streamline  increased  that  time  to  16  seconds.  Using  the  error  control 
scheme  with  a  ten  percent  tolerance  on  the  truncation  error  yielded  an 
integration  time  of  13.75  seconds  (for  two  streamlines)  while  a  one 
percent  tolerance  increased  the  integration  time  to  23  seconds. 
Integration  of  the  species  continuity  equations  after  each  time  step, 
with  two  streamlines  and  five  intermediate  points,  yielded  an  accep¬ 
table  solution  in  an  overall  computational  time  of  nearly  one  hour. 
However,  for  this  system,  it  is  possible  to  integrate  the  species 
continuity  equations  as  infrequently  as  once  every  30  time  planes. 


Figure  17.  Total  enthalpy  error  profiles  for  the  C-H-O-N  system. 


This  reduces  the  computational  time  to  approximately  five  minutes. 


C-H-O-N  With  C,^  System 

It  is  possible  that  unburned  fuel  may  be  present  in  the  reactive 
mixture  at  the  nozzle  inlet  due  to  incomplete  combustion  in  the  ramjet. 
A  five  percent  mass  fraction  of  propane  (CgHg)  was  added  to  the  C-H-O-N 
chemistry  system  described  above  (for  the  same  nozzle  geometry)  to 
investigate  the  effect  of  the  unburned  fuel  and  to  provide  initial 
conditions  at  the  nozzle  inlet  which  were  not  the  result  of  an  equi¬ 
librium  calculation.  All  mass  fractions  of  the  species  present  in 
the  C-H-O-N  system  were  reduced  by  5  percent  to  allow  for  the  incor¬ 
poration  of  propane.  The  solution  with  the  addition  of  propane  did 
not  differ  greatly  from  the  C-H-O-N  solution;  the  mass-averaged  tem¬ 
perature  at  the  throat  increased  by  approximately  70°R  while  the 
mass-averaged  velocity  magnitude  decreased  by  55  ft/sec.  The  mass 
fraction  of  propane  remained  almost  constant  from  the  nozzle  inlet 
to  the  exit.  The  addition  of  one  chemical  specie  and  one  reaction 
slightly  increased  the  computational  time  required  for  integration  of 
the  species  continuity  equations  (about  0.5  seconds  per  integration). 

H-F  System 

The  H-F  system  was  analyzed  because  it  contains  only  five  chemi¬ 
cal  species  and  six  reactions;  relatively  short  computational  times 
per  time  plane  can  be  achieved  with  this  system.  Also,  since  the 
C-H-O-N  system  results  were  near  the  frozen  solution,  the  H-F  system 
(at  a  high  temperature)  provided  the  opportunity  to  investigate  a 
system  near  equilibrium.  A  small  nozzle  geometry  was  chosen  for  this 


case  so  that  some  departure  from  equilibrium  might  occur  before  the 
throat.  As  shown  in  Appendix  E,  thp  nozzle  is  circular-arc  conical 
with  a  two  inch  inlet  radium,  a  one  inch  throat  radius,  a  45  degree 
angle  in  the  converging  section,  and  a  15  degree  angle  in  the  diver¬ 
gence.  Figure  18  shows  the  results  of  the  H-F  system  analysis.  The 
mass-averaged  temperature  profile  falls  between  the  equilibrium  and 
frozen  solutions  and  matches  the  ODK  [Reference  (10)]  solution  very 
well.  Note  that  this  system  is  near  equilibrium  upstream  of  the 
throat  but  departs  significantly  from  equilibrium  downstream  of  the 
throat.  A  difference  of  approximately  one  percent  in  the  inlet  and 
throat  mass  flows  occurred  for  both  the  isentropic  (y  =  constant)  and 

kinetics  solutions.  The  species  mass  fractions  deviated  by  less  than 
-10 

lx  10  along  the  center  and  wall  streamlines.  Figure  19  presents 
the  total  enthalpy  profiles  for  the  H-F  system.  For  this  particular 
problem,  the  solutions  for  the  fluid  dynamic  variables  along  the 
centerline  (upstream  of  the  throat)  were  still  oscillating  slightly 
after  more  than  700  time  planes  for  both  the  kinetics  and  isentropic 
solutions.  This  accounts  for  the  sign  difference  in  the  total 
enthalpy  error  along  the  center  in  Figure  19.  Figure  20  illustrates 
the  two-dimensional  character  of  the  solution  for  the  H-F  system. 

Note  the  significant  differences  in  temperatures  and  species  mass 
fractions  between  center  and  wall  values  which  develop  in  the  solu¬ 
tion.  Also  note  the  compression  along  the  wall  just  downstream  of 
the  throat.  This  occurs  at  the  junction  between  the  circular  arc  of 
the  throat  region  and  the  conical  divergence.  The  effect  is  present 
in  both  the  isentropic  (y  =  constant)  and  kinetics  solutions. 


Figure  18.  Temperature  profiles  for  the  H-F  system. 


enthalpy  error  profiles  for  the  H-F  system. 


For  the  H-F  system,  with  25  mesh  points  along  the  axis,  two 
chemical  kinetics  streamlines,  and  five  intermediate  points  specified, 
each  integration  of  the  species  continuity  equations  required  five 
seconds  on  the  CDC  6500  computer.  The  use  of  four  streamlines 
increased  that  time  to  10  seconds  per  integration.  Using  the  error 
control  scheme  with  a  ten  percent  tolerance  on  the  truncation  error 
yielded  an  integration  time  of  5.8  seconds  (two  streamlines)  while  a 
one  percent  tolerance  increased  the  integration  time  to  10.4  seconds. 
By  adjusting  the  time  step  multiplier  A  [see  equation  (51)]  it  was 
possible  to  integrate  the  species  continuity  equations  only  every- 
other  time  plane  for  the  H-F  system.  However,  after  more  than  800 
time  planes,  the  values  of  the  fluid  dynamic  variables  at  the  center- 
line  were  still  oscillating  with  little  evidence  of  damping.  Inte¬ 
grating  the  species  continuity  equations  after  each  time  step  yielded 
an  acceptable  solution  in  a  total  computational  time  of  one  hour. 

3.  RECOMMENDATIONS  CONCERNING  COMPUTATIONAL  TIME 

Recall  from  the  discussion  of  the  overall  numerical  algorithm 
that  the  solution  to  the  subject  problem  is  advanced  in  a  two-step 
procedure;  the  fluid  dynamic  variables  are  advanced,  and  then  the 
species  continuity  equations  are  integrated  along  streamlines  through 
the  flow.  The  second  step  is  much  more  time  consuming  than  the  first 
and  thus,  any  attempt  to  reduce  or  minimize  computational  time  should 
focus  on  the  integration  of  the  species  continuity  equations.  The 
most  significant  factors  regarding  the  time  required  for  this  inte- 


1)  the  number  of  chemical  species  and  reactions  in  the  chemical 
kinetics  model . 

2)  the  use  of  the  truncation  error  control  scheme  for  integration 
of  the  species  continuity  equations. 

3)  the  number  of  chemical  kinetics  streamlines  used. 

4)  the  frequency  of  integration  of  the  species  continuity  equa¬ 
tions  as  the  solution  is  marched  forward  in  time. 

Controls  for  these  four  factors  are  available  to  the  program  user  but 
experience  to  date  indicates  that  the  proper  set  of  controls  is  prob¬ 
lem  dependent.  Some  guidance  is  provided  in  the  following  discussion. 

The  chemical  kinetics  model  should  be  as  simple  as  possible. 

The  technique  of  screening  reactions  for  inclusion  in  the  chemical 
kinetics  model,  as  presented  in  Reference  (38),  is  recommended. 

The  error  control  scheme  is  valuable  for  setting  up  the  trunca¬ 
tion  error  control  parameter  NINT  and  for  refining  the  solution,  but 
it  is  generally  too  costly  to  use  throughout  an  entire  calculation. 

The  user  should  begin  the  analysis  by  computing  several  time  planes 
with  error  control  and  a  relatively  loose  tolerance  (5-10%).  Then, 
recompute  the  same  time  steps  with  a  fixed  number  of  intermediate 
points.  Adjust  that  fixed  number  of  points  until  satisfactory  agree¬ 
ment  with  the  solution  using  error  control  is  achieved.  The  error 
control  scheme  should  also  be  used  to  refine  the  solution  after  con¬ 
vergence  with  a  fixed  number  of  intermediate  points. 

For  the  cases  studied,  the  addition  of  chemical  kinetics  stream¬ 
lines  to  a  two  streamline  solution  does  not  make  dramatic  changes  in 


the  computed  results.  The  use  of  more  than  two  streamlines  signifi¬ 
cantly  increases  execution  time  and  thus,  should  be  reserved  for 
refining  the  solution. 

The  C-H-O-N  system  results  illustrated  the  dramatic  computational 
time  savings  which  can  be  realized  if  it  is  not  necessary  to  integrate 
the  species  continuity  equations  after  each  time  step.  Again,  some 
experimentation  in  this  regard  should  be  performed  for  each  new  problem. 


SECTION  VIII 


CONCLUSIONS 

A  numerical  method  has  been  developed  for  solving  the  equations 
governing  two-dimensional,  unsteady,  chemically  reacting  flow  in 
propulsive  nozzles.  The  unsteady  equations  are  hyperbolic  in  sub¬ 
sonic,  transonic,  and  supersonic  flow  regions  and  thus  can  be  solved 
by  well  developed  marching  techniques.  The  steady  state  solution  is 
obtained  as  the  asymptotic  solution  to  the  unsteady  equations,  with 
steady  flow  boundary  conditions  applied,  for  large  time.  The  overall 
numerical  algorithm  .  s  inconsistent  in  time  in  the  treatment  of  the 
species  continuity  equations  but  becomes  consistent  at  the  steady 
state  limit.  Interpolation  for  species  mass  fractions  is  not  required 
as  part  of  the  scheme  for  integration  of  the  species  continuity 
equations. 

The  primary  contribution  of  this  research  is  the  development  of  a 
production-type  computer  program  suitable  for  application  to  a  variety 
of  nozzle  problems.  Verification  of  the  results  of  the  subject  analy¬ 
sis  for  two-dimensional  subsonic  and  transonic  flows  is  difficult 
because,  at  present,  there  are  no  other  methods  available  for  the 
solution  of  this  problem.  However,  for  the  cases  studied,  the  solu¬ 
tions  are  quite  reasonable  and  agree  well  in  a  mass-ai.eraged  sense 
with  accepted  one-dimensional  finite-rate  chemical  kinetics  solutions. 


81 


Also,  when  applied  to  two-dimensional  supersonic  flow,  the  results 
of  the  subject  analysis  compare  favorably  with  an  accepted  method- 
of -characteristics  solution. 

Computational  time  has  been  found  to  be  highly  problem  dependent. 
It  is  especially  sensitive  to  the  number  of  chemical  species  and  reac¬ 
tions  in  the  problem  model,  the  use  of  the  truncation  error  control 
scheme  in  the  integration  of  the  species  continuity  equations,  the 
number  of  chemical  kinetics  streamlines  used,  and  the  frequency  with 
which  the  species  continuity  equations  are  integrated  through  the 
flow  field  as  the  solution  is  marched  forward  in  time.  Appropriate 
controls  for  these  factors  are  available  to  the  program  user.  For 
the  cases  which  have  been  analyzed  to  date,  solutions  have  been 
obtained  in  times  ranging  from  five  minutes  to  one  hour  using  a  CDC 
6500  computer. 

The  primary  value  of  t.e  subject  technique  is  that  it  provides 
the  analysis  of  the  subsonic  and  transonic  portions  of  the  nozzle  flow 
field,  including  two-dimensional  and  finite-rate  chemical  kinetics 
effects,  and  starting  from  nonuniform,  nonequilibrium  conditions  at 
the  nozzle  inlet.  Also,  it  provides  consistent  chemical  kinetics  data 
in  the  supersonic  flow  region  which  can  be  used  to  generate  the 
initial-value  line  (surface)  for  accurate,  multidimensional  method-of- 
characteristics  techniques.  The  method  of  the  subject  research,  when 
coupled  with  a  method-of-characteri sties  scheme  for  the  nozzle  diver¬ 
gence,  should  provide  a  highly  accurate  analysis  of  the  entire  nozzle 
flow  field. 
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APPENDIX  A 
GOVERNING  EQUATIONS 

1.  STATEMENT  OF  THE  GOVERNING  EQUATIONS 

The  governing  equations  for  fluid  flow  within  the  nozzle  of  a 
propulsion  system  are  derived  from  the  application  of  the  law  of  con¬ 
servation  of  mass,  Newton's  second  law  of  motion,  and  the  first  law 
of  thermodynamics  to  an  appropriate  fluid  model.  For  analyses  per¬ 
formed  for  the  purpose  of  performance  prediction,  the  assumption  of 
a  continuous,  inviscid,  nondiffusing,  and  adiabatic  fluid  is  generally 
adequate  to  accurately  model  the  flow. 

The  macroscopic  conservation  equations,  which  can  be  developed 
from  the  laws  stated  above,  are  derived  in  detail  in  many  basic  gas 
dynamics  texts  (Ref.  35  for  example).  These  equations  are  merely 
stated  here.  Applying  the  law  of  conservation  of  mass  to  the  flow 
through  a  fixed  volume  element  yields  the  continuity  equation. 

J£  +  V  •  (pV)  =  0  (A-l) 

where  p  is  the  density  of  the  fluid  and  V  is  the  velocity.  When 
Newton's  second  law  is  applied  to  an  element  of  unit  mass  moving  with 
the  fluid  (for  the  fluid  model  described  above),  Euler's  equations  of 
motion  are  derived. 


where  P  is  the  fluid  pressure,  and  for  unsteady,  two-dimensional 
flow 

D(  )  _  ILL  +  u3(  )  +  v  ii J. 

Dt  3t  3X  3y 

This  vector  equation  provides  two  scalar  component  equations  in  a 
two-dimensional  flow.  Note  that  body  forces  have  been  assumed  to  be 
negligible  in  equation  (A-2).  Finally,  the  energy  equation  is  derived 
from  the  application  of  the  first  law  of  thermodynamics  to  an  adia¬ 
batic  element  of  unit  mass  moving  with  the  fluid.  One  possible  form 
of  this  equation  is  as  follows  (Ref.  25). 

D(h  +  V2/2)  liP.Q  / «  o\ 

Dt  p  3 1  U 

where  h  is  the  fluid  enthalpy.  This  form  of  the  energy  equation 
shows  that  for  steady  flow,  the  total  enthalpy  (h  +  V  / 2)  is  constant 
along  streamlines.  This  fact  serves  as  a  useful  check  on  the  validity 
of  computed  results  for  steady  flows. 

In  a  two-dimensional  flow,  equations  (A-l)  to  (A-3)  provide  four 
scalar  equations  In  the  five  unknowns  u,  v,  p,  p,  and  h,  where  u  and 
v  are  the  velocity  components.  In  equilibrium  and  frozen  flows,  an 
equation  of  state  of  the  form  h  *  h(p,p)  is  added  to  provide  a  fifth 
equation  and  to  complete  the  formulation  of  the  problem.  In  flows 
with  chemical  nonequilibrium,  equations  (A-l)  to  (A-3)  are  still  valid 
provided  that  the  thermodynamic  variables  are  given  an  "extended" 
definition  (Ref .  25).  It  is  assumed  that  the  nonequilibrium  system 
of  fixed  volume  is  in  mechanical  and  thermal  equilibrium  and  therefore 
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has  a  definite  pressure  and  temperature.  Also,  it  is  assumed  that 
the  system  is  homogeneous  in  space  and  therefore  its  composition 
can  be  specified  by  giving  the  number  of  moles  of  each  chemical 
species  present.  Finally,  it  is  assumed  that  pressure,  temperature, 
and  the  number  of  moles  of  each  of  the  chemical  species  are  not 
independent,  but  are  related  to  one  another  through  an  equation  of 
state.  To  insure  that  the  equations  of  state  for  a  nonequilibrium 
situation  reduce  to  the  equilibrium  form  when  the  equilibrium  compo¬ 
sition  is  substituted  into  those  equations,  the  following  requirement 
is  made  (Ref.  25).  "The  state  equation  for  any  property  of  a  system 
in  chemical  nonequilibrium  as  a  function  of  any  other  two  properties 
and  all  the  mole  fractions  of  the  constituent  species  is  identical  in 
form  to  the  corresponding  equation  for  the  system  in  thermodynamic 
equilibrium."  For  example,  it  is  this  requirement  which  allows  the 
use  of  the  familiar  thermal  equation  of  state  for  a  mixture  of  perfect 
gases  in  a  flow  with  chemical  nonequilibrii 

Consider  an  equation  of  state  of  the  following  form  for  the 
nonequilibrium  flow. 

h  =  h(p,P,C1,...,Cn) 

where  (i  =  l . n)  are  the  mass  fractions  of  each  of  the  n  chemi¬ 

cal  species  in  the  reacting  fluid.  The  addition  of  this  equation  of 
state  to  equations  (A-l)  to  (A-3)  provides  a  set  of  five  equations 
in  five  plus  n  unknowns.  All  additional  n  equations  are  provided  by 
the  application  of  the  law  of  conservation  of  mass  to  the  flow  of 
each  of  the  chemical  species  through  a  fixed  volume  element.  This 
yields  the  species  continuity  equations. 
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where  p.  is  the  mass  density  of  species  i  and  is  the  species 
source  function. 

As  suggested  above,  it  will  be  assumed  in  this  analysis  that  the 
working  fluid  is  a  mixture  of  thermally  perfect  gases.  Then  the 
thermal  equation  of  state  is: 

n 

P  *  pT  Z  C.R.  =  pRT  (A-6) 

i=l  1  1 


The  caloric  equation  of  state  is: 


n  T 

h  *  Z  C.h.  where  h.  =  /  C  <  dT  +  h.°  (A-7) 

4-1  i  i  i  T  Pi  i 

o 


In  equations  (A-6)  and  (A-7),  R^  are  the  species  gas  constants,  R  is 
the  gas  constant  for  the  mixture,  h^  are  the  species  enthalpies, 

Cpi  are  the  species  specific  heats  at  constant  pressure,  TQ  is  the 
reference  temperature,  T  is  the  fluid  temperature,  and  h^°  are  the 
species  energies  of  formation. 


Summary  of  major  assumptions  and  governing  equations 
Major  assumptions: 

-  continuum 

-  inviscid 

-  adiabatic 

-  nondiffusing 

-  negligible  body  forces 

-  chemical  nonequilibrium  only 

-  fluid  is  a  mixture  of  thermally  perfect  gases 


Governing  equations:  (vector  form) 


Global  continuity: 

+  V  •  (pV)  =  0 

(A-l) 

Momentum: 

o 

ii 

Q. 

> 

+ 

l>|+J 

0|Q 

CL 

(A-2) 

Energy: 

D(h  +  V2/2)  1  3  P _ 

Dt  P9t  u 

(A-3) 

Species 

continuity: 

+  V  •  (p.v)  =  oi  (i =  l,...,n) 

(A-5) 

Thermal  equation 

n 

(A-6) 

of  state: 

P  -  PT  E  C.R. 
i=l  1  1 

Caloric  equation  n  T 

of  state:  h  =  Z  C.h.  where  h.  =  /  C  .  dT  +  h.°  (A-7) 

l'=1  '  '  1  T  P*  1 
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2.  REARRANGEMENT  OF  THE  SPECIES  CONTINUITY  EQUATIONS  AND  THE  ENERGY 
EQUATION 

The  species  continuity  equations  (A-5)  can  be  simplified  by 
expressing  them  in  terms  of  species  mass  fractions  rather  than  species 
mass  densities.  Substituting  C^  =  p^/p  into  equation  (A-5)  and  using 
the  global  continuity  equation  (A-l)  yields: 

DC  j  o1 

=  —  (i  =  l,...,n)  (A-8) 

In  this  characteristic  form,  the  species  continuity  equations  are 
one-dimensional  equations  which  apply  along  a  particle  path.  They 
are  integrated  in  this  form  in  the  computational  scheme. 


The  energy  equation  (A-3)  can  also  be  manipulated  into  a  form 
which  is  more  convenient  for  analysis.  In  this  form  the  species 
source  function  appears  directly.  The  species  source  function 
contains  within  it  the  set  of  chemical  reactions  and  appropriate 
reaction  rates  for  the  assumed  chemistry  model  (see  Appendix  B). 

By  forming  the  dot  product  of  the  velocity  vector  V  with  equation 
(A-2),  it  is  possible  to  develop  an  expression  for  . 


pv  *  +  v  ‘  VP 


-  n  (¥-)  +  21  IP 

-  pD(  2  )  +  Dt  -  3t 


(A-9) 


2 

Solving  equation  (A-9)  for  D(V  /2)  and  substituting  into  equation 
(A-3)  yields 


Dh  1  =  n 

Dt  '  p  Dt  u 


(A- 10) 


Introducing  the  thermal  and  caloric  equations  of  state,  equations 
(A-6)  and  (A-7),  respectively,  and  the  species  continuity  equation 
(A-8),  yields  the  desired  form  of  the  energy  equation.  Taking  the 
substantial  derivative  of  equation  (A-7)  gives 


Dh 


n 


Dt  ~  ^ 
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(A-ll) 

Logarithmic  differentiation  of  equation  (A-6)  provides  an  expression 
for  DT/Dt  in  equation  (A-ll).  Equation  (A-ll)  can  then  be  used  to 
replace  Dh/Dt  in  equation  (A-10).  The  result  is 
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(A- 12 ) 


OP 

Dt 


Dfi 

Dt 


where 


*k  =  tYfRiT  "  Wf-Dhiloi 


(A- 13 ) 


and  and  y^  are  the  frozen  speed  of  sound  and  ratio  of  frozen  speci¬ 
fic  heats,  respectively. 


3.  THE  GOVERNING  EQUATIONS  FOR  TWO-DIMENSIONAL  PLANAR  OR  AXISYMMETF 
FLOW 

The  geometries  of  interest  in  this  research  include  two-dimen 
sional  planar  and  axisymmetric  flows.  The  equations  governing  two- 
dimensional  axi symmetric  flow  are  normally  expressed  in  cylindrical 
coordinates,  but  they  can  be  transformed  to  the  notation  of  a  Cartesian 
coordinate  system  by  means  of  the  following  transformation  equations. 


r  =  y 

z  =  x 


Comparing  the  equations  for  two-dimensional  planar  flow  with  the 
transformed  set  of  equations  [by  equations  ( A-14) ]  for  a  two-dimensional 
axisymmetric  flow,  it  is  seen  that  the  equations  are  identical  except 
for  the  additional  term  pv/y  in  the  global  continuity  equation  for 
axisymmetric  flow.  Therefore,  the  following  set  of  equations,  in 
Cartesian  coordinate  notation,  applies  to  both  two-dimensional  planar 
and  two-dimensional  axisynmetric  flows. 

Pt  +  Upx  +  vpy  +  pux  +  pVy  +  =  0  (A- 15) 
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ut  +  uux 

+ 

vuy 

♦  px/p  = 

0 

(A- 16 ) 

vt  +  uvx 

+ 

Wy 

*  v° = 

0 

(A-17) 

Pt  +  UPX 

+ 

vPy 

-  af2(pt 

+  upx  +  vpy>  *  *k 

(A-18) 

a. 

Cit  +  uC 

ix 

+  vC 

<<’ 

ii 

^1- 

(i  =  1, . . .  ,n) 

(A-19) 

where  the  subscripts  denote  partial  differentiation.  In  equation 
(A-15),  e  is  zero  for  planar  flow  and  one  for  axisymmetric  flow. 

4.  TRANSFORMATION  OF  THE  GOVERNING  EQUATIONS  TO  THE  COMPUTATIONAL  PLANE 

Since  the  interior  points  in  this  analysis  are  to  be  treated  by 
a  fixed  grid  technique,  it  is  convenient  to  transform  the  physical 
(x,y,t)  plane  to  a  rectangular  (5,0,1)  plane  in  which  the  differencing 
is  performed.  The  following  coordinate  transformation  is  used  (see 
Figure  A-l) 

T  =  t 

(A-20) 

C  =  x 

y  -  yc(x) 

n'  -  *c<*> 

where  y  (x )  and  yul(x)  represent  the  nozzle  centerbody  and  wall  coor- 

L  W 

di nates,  respectively,  as  functions  of  x.  Application  of  this  trans¬ 
formation  yields: 

ill.  9JJ. 

3t  "  3x 

1LA  =  LLL  +  a  ILL 

3x  85  3  0 

ill  =  B  3(  ) 


( A- 21) 
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an.  6^c  n6(a^w  ^ 

-8  — -^(^x  3X) 


g  =aa 


ay  "  yw(x)  -  yc(x) 


For  convenience,  the  following  two  variables  are  defined: 


v  i  au  +  8v  (A-22) 

n  =  yc  +  n/8  (A-23) 

Applying  the  transformation  in  (A-21)  to  equations  (A-15)  to  (A-19) 
and  introducing  the  definitions  given  by  equations  (A-22)  and  (A-23) 
yields  the  following  forms  of  the  continuity,  momentum,  energy,  and 
species  continuity  equations,  which  are  used  in  this  analysis. 

Summary  of  the  governing  equations  iji  the  computational  plane 


pt  +  upt,  +  ^Pn  +  puc  +  potun+p^vn  +  ePv/^  =  0 

(A-24) 

uT  +  uu^  +  vun  +  P^/p  +  ap^/p  =  0 

(A-25) 

vT  +  uv^  +  Vvn  +  ep^/p  =  0 

(A- 26) 

PT  +  uP?  +  vPn  -  af?(PT  +  upc  +  vPn)  = 

(A— 27 ) 

ciT  +  ueic  +  vCin  =  ai/P  (i  =  l, . . .  ,n) 

(A-28) 

n 

P  =  pT  Z  C.R.  =  pRT 
i=l  1  1 

(A-6) 

n  T 

h  =  E  C.h.  where  h.  =  /  C_.dT  +  h.° 
i=l  11  To  P 

(A-7) 
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Note  that  equation  (A-28)  may  be  written  as 


(i =  l,...,n) 


where  DC^/Dx  is  the  change  in  following  a  particle  path  in  the 
computational  plane. 
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APPENDIX  B 


THE  CHEMICAL  KINETICS  MODEL 

1.  REACTION  MECHANISM  AND  THE  SPECIES  SOURCE  FUNCTION 

In  Appendix  A,  it  was  noted  that  for  the  nonequilibrium  flow  of 
a  gas  mixture,  a  third  variable  {the  species  concentrations)  must  be 
introduced  into  the  equation  of  state  to  specify  the  chemical  state  of 
the  reactive  mixture.  This,  in  turn,  required  the  addition  of  n  species 
continuity  equations  to  complete  the  set  of  governing  equations.  In 
this  section,  an  equation  for  the  source  term  in  the  species  continuity 
equations  is  presented. 

Implicit  within  the  equation  for  the  species  source  function  is 
a  reaction  mechanism  (i.e.,  a  set  of  coupled  elementary  reactions  whose 
reaction  rates  are  functions  of  temperature  only).  As  an  example  of 
a  reaction  mechanism,  consider  the  overall  reaction  between  hydrogen 
and  flourine: 

H2  +  «j>F2  =  aH2  +  bF2  +  cHF  +  dH  +  eF  (B-l) 

where  4>  is  the  molar  oxidizer  to  fuel  ratio,  and  a  through  e  are  the 
moles  of  the  various  species  in  the  completed  reaction.  The  science 
of  chemical  kinetics  is  concerned  with  the  determination  of  a  suitable 
reaction  mechanism  for  overall  reactions  like  that  specified  by  equa¬ 
tion  (B-l).  The  reaction  mechanism  should  give  an  overall  effect 


which  is  in  agreement  with  experimental  observations.  For  the  overall 
reaction  given  by  equation  (B-l),  Cherry  (36)  has  proposed  the  follow¬ 
ing  reaction  mechanism: 


F2  +  M  t  IF  +  M  (B-2a) 
H2  +  M  %  2H  +  M  (B-2b) 
HF+M^H  +  F  +  M  (B-2c) 
HF  +  F  t  H  +  F2  (B-3a) 
HF  +  H  X  H2  +  F  (B-3b) 
2HF  t  H2  +  F2  (B-3c) 


Reactions  (B-2)  are  called  dissociation-recombination  (or  third  body) 
reactions,  while  reactions  (B-3)  are  known  as  binary  exchange  reac¬ 
tions.  The  symbol  M  in  equation  (B-2)  represents  the  third  body 
involved  in  the  reactive  collision  and  it  can  be  any  of  the  chemical 
species  present  in  the  reaction  mechanism.  A  general  reaction  equa¬ 
tion  which  will  be  used  to  represent  any  reaction  mechanism  is 


n 

£ 

i=l 


n 

£ 

i=l 


^  xi 


(  J  ““  1  9  *  •  • 


(B-4) 


i  ii 

where  v^j  and  are  the  stoichiometric  coefficients  of  the  reactants 

and  products  respectively,  denotes  the  ith  chemical  species,  and 

k-j  and  k„.  are  the  forward  and  reverse  reaction  rate  coefficients, 
tj  ro 

respectively,  for  the  jth  reaction  of  the  m  reactions  in  the  reaction 
mechanism. 


The  law  of  mass  action  (Ref.  25)  states  that  "the  rate  at  which 
an  elementary  reaction  proceeds  is  proportional  to  the  product  of  the 
molar  concentrations  of  the  reactants  each  raised  to  a  power  equal  to 
its  stoichiometric  coefficient  in  the  reaction  equation."  Using  this 
law  together  with  equation  (B-4),  it  is  possible  to  develop  an  expres¬ 
sion  for  the  net  production  of  the  ith  chemical  species  due  to  the  jth 
reaction  in  the  reaction  mechanism  (Ref.  31  ).  Thus, 


d[X.] 

-a r  J 


n  vn 
kr1  n  [X,]  ,J] 
rj  i=1  i 


(B-5) 


Converting  equation  (B-5)  to  a  mass  basis  and  summing  over  all  m  reac¬ 
tions  in  the  reaction  mechanism  yields  the  equation  for  the  species 
source  function. 


a. 
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_  m  „ 
m.  £  (v 
1  j=l  1J 


1  n  or  vi  i 

V  [kfj  "  (Pri>  3  -  krj 

TJ  i=l  m ^  rj 


n  pC.  v. 
n  (-^1)  1J] 

i=l  V 


(B-6) 


where  m^  is  the  molecular  weight  of  the  ith  species.  Note  that  the 
species  source  function  is  the  forcing  term  in  the  species  continuity 
equations  (A-8)  and  that  a.  is  also  present  in  the  forcing  term  of  the 
energy  equation  (A-13).  Therefore,  equation  (B-6)  provides  the  link 
between  the  chemical  kinetics  model  and  the  governing  equations  for 
the  problem  of  interest. 


2.  REACTION  RATE  COEFFICIENTS 


I 


Reaction  rate  coefficients  are  not  readily  predicted  at  the 
present  time  and  therefore,  they  are  generally  found  experimentally. 
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The  form  for  the  reverse  reaction  rate  coefficient  which  is  used  in 


this  analysis  is 

krj  =  aij  T  exP('bj/]rT)  (B-7) 

where  a..,  n.,  and  b.  are  empirical  coefficients,  R  is  the  universal 
gas  constant,  and  T  is  the  local  gas  temperature.  The  forward  and 
reverse  reaction  rate  coefficients  for  a  given  reaction  are  not  indepen¬ 
dent  but  are  coupled  through  the  equilibrium  constant  for  that  reaction. 

Using  the  equilibrium  constant  based  on  partial  pressures  (K  ,),  it 

P>  J 

is  found  that  (Ref.  31) 

kf  j  _  v  -  Av . 

krj  "  Kp,j  (AT)  J  (B-8) 

where 

n  ii  i 

Av.  =  l  (v,  -  v.)j  (B-9) 

J  1=1  1  1 

Note  that  the  difficulties  in  making  accurate  measurements  for  the 
determination  of  reaction  rate  coefficients,  the  necessity  of  extrapo¬ 
lating  experimental  data  from  one  situation  to  another,  and  the  uncer¬ 
tainty  as  to  the  reaction  mechanism,  are  all  significant  factors 
limiting  the  accuracy  of  the  analysis  presented  in  this  research.  How¬ 
ever,  experience  with  one-dimensional  analyses  of  reacting  flows 
indicates  that  the  chemical  kinetics  model  used  here  is  adequate  for 
performance  prediction. 
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3.  TREATMENT  OF  HYDROCARBONS 


Unburned  hydrocarbons  may  be  present  in  the  gas  mixture  at  the 
combustor  exit  and  therefore,  provision  for  these  species  must  be  made 
in  the  chemical  kinetics  model.  Edelman,  et.  al.  (26)  and  Edelman  and 
Harsha  (27)  have  developed  a  kinetics  model  for  use  in  combustor 
analyses  which  includes  a  "sub-global"  partial  oxidation  step: 

CnHm  +  t  °2  *  I  H2  +  nCO  (B-10) 

Note  that  this  reaction  proceeds  in  the  forward  direction  only.  The 
reaction  rate  coefficient  for  equation  (B-10)  is  determined  empirically 
and  is  given  by  (Ref.  27): 

a.  =  6.0  x  1014  P0-3  cj  h  cn  (B-ll) 

3  Vm  °2 

where  P  has  units  of  atmospheres  and  c  represents  species  concentra¬ 
tions  in  g-moles/cc. 

4.  THE  SPECIES  SOURCE  FUNCTION  FOR  DISSOCIATION-RECOMBINATION  REAC¬ 
TIONS  AND  BINARY  EXCHANGE  REACTIONS 

Equation  (B-6)  for  the  species  source  function  is  simplified  for 
dissociation-recombination  and  binary  exchange  reactions  in  the  follow¬ 
ing  discussion. 

Dissociation-Recombination  Reactions 

Dissociation  recombination  reactions  are  of  the  form 

A  +  M^B  +  C  +  M  (B-12) 

From  equation  (B-9),  Av.  =  1  for  all  these  reactions.  Substituting 
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equation  (B-8)  into  equation  (B-6)  yields 
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(B-13) 


(B-14) 


(B-15) 


The  effect  of  different  third  bodies  is  taken  into  account  through  the 
use  of  third  body  reaction  rate  ratios  M.  (see  section  4). 

Binary  Exchange  Reactions 

Binary  exchange  reactions  are  of  the  form 

A  +  BJC  +  0  (B-16) 

In  this  case,  Av.  =  0.  Thus,  equation  (B-6)  becomes 

J 
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(B-18) 


where 


K  •  =  Kn  , 


n  - 

n  (m.)  1J 

7-1  1 


J  P  »  J  Q  _  Vii 

B-E  2  (m. )VlJ 

i=l  1 


(B-19) 


The  Hydrocarbon  Reaction 

The  hydrocarbon  oxidation  as  shown  in  equation  (B-10)  can  be 
treated  as  a  special  case  of  a  binary  exchange  reaction.  Rewriting 
equation  (B-10)  in  the  form  of  (B-16),  and  so  that  it  proceeds  from 
right  to  left,  gives 

0 

?  H2  +  n  CO  t  CnHm  +  \  02  (B-20) 

Since  K„  .  is  equal  to  the  ratio  of  the  forward  and  reverse  reaction 
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rates,  then  K,,  .  is  zero  for  the  reaction  in  equation  (B-20).  There- 

■  »J 

fore,  the  species  source  function  for  the  special  case  of  equation 
(B- 10)  becomes 

_  o  n  .  n  k  . 

o  =  m  p2  (v  -  v  )  [-  n  (C.)  1J]  — - -  (B-21) 

11  *J  'J  -isi  1  n  _ 

11  n  (m.)  13 

i=l  1 


5.  THIRD  BODY  REACTION  RATE  RATIOS 

Each  dissociation-recombination  reaction  has  a  different  reverse 
reaction  rate  depending  upon  the  particular  third  body  involved  in  the 
collision.  This  effect  is  taken  into  account  in  the  chemical  kinetics 
model  by  the  use  of  third  body  reaction  rate  ratios.  It  is  assumed 
that  the  overall  reverse  reaction  rate  coefficient  for  the  jth  disso¬ 
ciation-recombination  reaction  is  given  by 


M .  k  .  =  E 


J  rj 


A=1 


C*  krtj 


(B-22) 


where  the  summation  over  g,  represents  all  possible  third  bodies  and 
is  the  mass  fraction  of  the  tth  third  body.  Implicit  in  this 
equation  is  the  assumption  that  the  overall  reverse  reaction  rate  is 
a  mass  weighted  average  of  the  reverse  reaction  rates  for  each  differ¬ 
ent  third  body.  According  to  References  (36)  and  (37),  the  temperature 
dependence  of  recombination  rates  can  be  reasonably  assumed  to  be 
independent  of  the  third  body  so  the  recombination  rate  associated 
with  the  tth  species  (third  body)  can  be  expressed  as 

kr*j  ‘  -d  T'"J  «P<-V*T)  (B-23) 
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where  is  dependent  on  the  third  body.  Substituting  this  expression 
into  equation  (B-18)  yields 


VrJ  '  J,  C*  T  J 


(B-24) 


£ 

or,  using  equation  (B-7), 


Therefore,  M^  is  defined  by 

M.  =  I  (^i)  C£  (B-27) 

J  £-1  aij 

Physically,  M.  is  the  mass  weighted  ratio  of  different  third  body 

J 

reaction  rate  premultipliers,  a  .,  to  the  known  reaction  rate  pre- 

X*  J 

multiplier  a... 

'  J 

6.  REACTION  MECHANISM  USED  IN  THIS  RESEARCH 

The  reaction  mechanism  used  in  this  research  is  the  same  as  that 
used  in  Reference  (3)  with  the  addition  of  the  hydrocarbon  reaction  dis¬ 
cussed  in  Section  3.  For  a  nonequilibrium,  chemically  reacting  flow 
of  a  system  of  thermally  perfect  gases  composed  of  the  six  elements 
carbon,  hydrogen,  oxygen,  nitrogen,  flourine,  and  chlorine.  Gold  and 
Weekly  (28)  have  shown  that  19  species  and  48  chemical  reactions  should 
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(B-25) 
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be  considered.  The  chemical  species  are  shown  in  Table  B-I  and  the 
48  reactions  are  shown  in  Table  B-II.  Note  that  13  reactions  are 
shown  in  Table  B-II.  Note  that  13  reactions  are  of  the  dissociation- 
recombination  type  while  the  remaining  35  are  binary  exchange  reac¬ 
tions. 

Experience  with  this  and  other  reaction  mechanisms  has  indicated 
there  may  be  cases  when  not  all  of  the  reactions  corresponding  to  a 
given  set  of  chemical  species  should  be  included  in  the  analysis.  The 
number  of  reactions  in  the  reaction  mechanism  is  a  significant  para¬ 
meter  affecting  execution  time.  Therefore,  no  more  reactions  than  are 
necessary  should  be  included  in  the  analysis.  Reference  (38)  provides 
a  grading  of  the  various  reactions  appropriate  for  CHON  and  HF  pro¬ 
pellant  systems  as  follows: 

A  -  Important  reaction  between  important  species  with  reliable 
rate  data  from  a  competent  review 

B  -  Energetically  less  significant  reaction  or  rate  data  of 
doubtful  quality 

C  -  Reaction  between  minor  species  which  may  be  of  significance 
in  some  nozzle  flows 

X  -  Possibly  significant  reaction  with  estimated  rate  data 
(not  recommended  for  use) 

Guidelines  for  the  selection  of  reactions  from  the  overall  graded  set 
of  reactions  are  provided  in  Reference  (38). 


TABLE  B-l 

CHEMICAL  SPECIES  CONSIDERED 


TABLE  B-2 


Chemical  Reaction 


Chemical  Reaction 


C  H  +  ^  0,  -  H,  +  nCO 
n  m  L  L  L  L 

HC1  +  HC1  t  H2  +  Cl2 

C02  +  M  X  CO  +  0  +  M 

HC1  +  0  X  OH  +  Cl 

H20  +  M  t  OH  +  H  +  M 

HF  +  Cl  X  HC1  +  F 

CO  +  MjC  +  O  +  M 

HF  +  F  X  H  +  F2 

Cl2  +  M  t  2C1  +  M 

HF  +  H  X  H2  +  F 

F2  +  M  t  2F  +  M 

HF  +  HF  X  H2  +  F2 

HC1  +  M  t  H  +  Cl  +  M 

HF  +  0  t  OH  +  F 

HF  +  MiH  +  F  +  M 

HF  +  OH  t  H20  +  F 

H2  +  M  t  2H  +  M 

H2  +  Cl  t  HC1  +  H 

N2  +  M  J  2N  +  M 

H2  +  0  t  OH  +  H 

NO  +  MtN+O  +  M 

H2  +  02  t  20H 

OH  +  MtO  +  H  +  M 

N?  +  0  t  NO  +  N 

02  +  M  t  20  +  M 

\D  f 

+ 

o 

+  + 

ro 

o 

C1F  +  M  t  Cl  +  F  +  M 

NO  +  H  t  N  +  OH 

C02  +  H  t  CO  +  OH 

NO  +  0  Z  N  +  02 

co2  +  0  t  CO  +  o2 

02  +  H  t  OH  +  0 

H20  +  Cl  t  OH  +  HC1 

Cl  +  C1F  t  Cl2  +  F 

H20  +  H  t  OH  +  H2 

F  +  C1F  t  Cl  +  F2 

H20  +  0  t  20 H 

HF  +  Cl  t  Cl F  +  H 

CO  +  CO  J  co2  +  C 

HC1  +  F  t  C1F  +  H 

CO  +  H  t  C  +  OH 

HC1  +  HF  t  C1F  +  H2 

CO  +  N  t  C  +  NO 

HF  +  C1F  t  F2  +  HC1 

CO  +  NO  t  CO,  +  N 

HF  +  Cl,  t  C1F  +  HC1 

CO  +  0  t  C  +  0, 

Cl F  +  Cl F  t  F,  +  Cl, 

HC1  +  Cl 


u 


APPENDIX  C 

REFERENCE  PLANE  CHARACTERISTIC  RELATIONS 

1.  GOVERNING  EQUATIONS 

The  equations  that  are  employed  in  this  analysis  to  model  unsteady, 
two-dimensional,  chemically  reacting,  nonequilibrium  flow  are  developed 
and  transformed  to  the  computational  U.n.x)  plane  in  Appendix  A. 


These  equations  are  restated  here  for  convenience. 

Pr  +  +  vpn  +  pu^  +  pau^  +  pSv^  +  epv/n  =0  (C-l) 

u  +  uu  +  vu  +  P  /p  +  cxP  / p  =  0  (C-2) 

t  c,  n  s»  n 

vx  +  uvc  +  ^Vn  +  ePr/p  =  0  (C-3) 

PT  +  uP^  +  vP^  -  af2(pT  +  UPC  +  7pn)  =  (C-4) 

Cit  +  uCic  +  ^Cin  =  °i/p  (i*l,...,n)  (C-5) 

h  =  h(P,p,C.)  (C-6) 


The  subscripts  x,  and  n  in  these  equations  denote  partial  differen¬ 
tiation  and  i  denotes  the  ith  chemical  species.  Equation  (C-6)  is  the 
general  form  of  the  equation  of  state  and  since  it  is  an  algebraic 
equation,  it  is  not  stated  again  in  the  remainder  of  this  appendix. 
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2.  CONSTANT  n  REFERENCE  PLANE 

When  the  method  of  characteristics  is  applied  to  unsteady,  two- 
dimensional  flow,  the  characteristic  surfaces  and  corresponding  com¬ 
patibility  equations  are  determined.  Hoffman  (39)  shows  that  the 
characteristic  surfaces  for  this  case  are  of  two  types;  stream  surfaces 
and  wave  surfaces.  The  envelope  of  the  stream  surfaces  is  the  pseudo¬ 
pathline  (the  particle  trajectory  in  W  space),  and  the  envelope  of 
the  wave  surfaces  is  the  Mach  cone.  The  curves  of  tangency  between  the 
wave  surfaces  and  the  Mach  cone  are  the  bicharacteristics.  One  com¬ 
patibility  equation  is  valid  on  each  wave  surface,  one  compatibility 
equation  is  valid  on  each  stream  surface,  and  one  compatibility  equa¬ 
tion  is  valid  along  the  pseudo-pathline.  Figure  C-l  illustrates  the 
Mach  cone,  bicharacteristics,  and  pseudo-pathline  as  discussed  above. 

A  reference-plane  characteristic  scheme  is  used  in  this  analysis  rather 
than  a  bicharacteristic  scheme  because  the  increased  complexity  and 
computational  time  of  the  bicharacteristic  scheme  is  not  warranted  in 
view  of  the  accuracy  limitations  associated  with  the  chemical  kinetics 
model  (see  Appendix  B). 

In  reference-plane  characteristic  schemes,  derivatives  with  res¬ 
pect  to  one  of  the  independent  variables  are  approximated  and  treated 
as  forcing  terms,  thus  reducing  the  number  of  independent  variables 
in  the  problem  by  one.  For  example,  in  the  constant  n  reference-plane 
scheme,  all  derivatives  with  respect  to  n  are  approximated  (by 
MacCormack's  method  or  some  other  suitable  method)  and  are  placed  on 
the  right-hand  side  of  the  equal  sign  in  equations  (C-l)  to  (C-5). 


Figure  C-l.  Mach  cone  and  pseudo-particle  path. 


f 


% 

i 


This  yields 

pt  +  up^+  pu^  =  Ijjj 

uT  +  uuc  +  P^/p  =  i|>2 
vT  +  uv^  =  ip3 

PT  +  Uf>C  ‘  af2(pT  +  Up^  =  ^4 

CiT  +  uC.^  =  (i=l,..  .,n) 

where 

^1  =  “vpn  "  potUn  "  p0vn  "  epV/^n 
*2  *  -*%  -  “’’n  /p 

*3  =  -“vn  '  SVP 

*4  *  -Vpn  *  af2*  %  +  ♦k 
*5  -  -Vcin  +  pi/p 

Characteristic  Curves 

A  linear  combination  of  the  governing  equations  can  be  formed  by 
multiplying  equations  (C-7)  to  (C-ll)  by  2,.  (j=l,2,...,4+n),  respec- 

J 

tively,  and  then  sumning  them.  The  notation  is  simplified  if  through 
*4+n  are  rePlacecl  ^  *-5*  The  resulting  linear  combination  is 

tj(px  +  up^  +  pu^  -  +  £2(ut  +  uu^  +  P^/p  -  <p2) 

+p,3(vt  +  uv^  -  i^3)  +  t4[Px  +  uP^  -  af2(px  +  up^)  -  i|>4]  (C-12) 

+*5(Cit  +  uC.c  -  *5)  -  0 


(C-7) 

(C-8) 

(C-9) 

(C-10) 

(C-ll) 
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Rearranging  equation  (C-12)  yields 

2  2 
pc(u*l  -  af  ut4)  +0^!  -  af  S.4) 

+  U^pJlj  +  uj^)  +  uT(  &2) 

+  v^(uJl3)  +  vt(ji3)  (C-13) 

+  P?U2/p  +  U£4)  +  Pt(A4) 

+  C^CuJlg)  +  C^Cig)  =  2^2  +  ^3^3  +  ^4  +  ^-5^5 

The  following  set  of  vectors  can  be  defined  where  the  components  are 
the  coefficients  of  the  partial  derivatives  in  equation  (C-13). 


V 

(u 

-  af2u£4  ,  -  af2£4) 

(C-14) 

V 

(p*! 

U&2  »  &  2) 

(C-15) 

V 

(UA3 

»  £3^ 

(C— 16) 

V 

,0. 

CVJ 

+  Ut4  ,  £4) 

(C— 17 ) 

V 

(u*5 

»  *5) 

(C-18) 

Now,  introduce  the  notation  that  dyiP  is  the  derivative  of  p  in  the 
direction  of  Wj,  d^u  is  the  derivative  of  u  in  the  direction  of  W,,, 
etc.  Recall  that  the  directional  derivative  of  a  function  f  in  the 
direction  of  W  is  found  by  taking  the  dot  product  of  the  gradient  of 
f  with  the  vector  W.  Therefore,  with  the  notation  above,  and  the 
vectors  defined  by  equations  (C-14)  to  (C-18),  equation  (C-13)  can  be 
written  as 


dWlp  +  dW2u  +  dW3v  +  dW4P  +  dW5  Ci  = 


(C-19) 


+  z2^2  +  ^a  +  £4^4  +  Z5^5 


If  the  (j=l,2,...,5)  can  be  chosen  so  that  the  vectors  W. 

J  J 

(j=l,2, . . . ,5)  are  linearly  dependent  (lie  in  one  direction),  then  the 
curve  which  contains  the  vectors  W.  is  called  the  "characteristic 

J 

curve,"  its  normal  is  called  the  "characteristic  normal,"  and  equa¬ 
tion  (C-19)  is  called  the  "compatibility  equation."  If  N  =  (N  ,N  ) 

C  T 

is  the  characteristic  normal  in  the  (c,r)  plane,  then  N  and  W.  must 

J 

be  related  by 


IT  •  V.  =  0  (j-1,2 . 5) 


(C-20) 


Expanding  equation  (C-20)  yields 


(u*j  -  af2u*4)Nc  +  (£j  -  af2£4)NT  =  0 


(C-21) 


(p*j  +  u«-2)N^  +  £2Nt  =  0 


(C-22) 


(u*3)Nc  +  £3Nt  -  0 


(C-23) 


(i2/p  +  U£4)NC  +  £4Nt  =  0 


(C-24) 


(u^5)nc  +  h5nt  =  0 


(C-25) 


Note  that  the  subscripts  c  and  x  in  equations  (C-21)  to  (C-25)  denote 
the  components  of  the  characteristic  normal  in  directions  of  z,  and  x, 
respectively.  In  matrix  form,  equations  (C-21)  to  (C-25)  can  be 
written  as  follows: 


UVNX 


uVNt 


0  -a‘(uN  +N  ^  0 


uVNt 


uNc+Nt 


unc+nt 


(C-26) 


Equation  (C-26)  is  a  system  of  homogeneous  equations  and  therefore,  if 
it  is  to  have  a  nontrivial  solution,  the  determinant  of  the  coeffi¬ 
cient  matrix  must  be  zero.  Setting  that  determinant  to  zero  yields 


(uN^  +  Nt)3[(uNc  +  Nt)2  -  af2Nc2]  =  0 


(C-27) 


Equation  (C-27)  has  two  possible  solutions.  Setting  the  first  factor 
to  zero  yields 


uN^  +  Nx  =  0 


(C-28) 


If  the  term  in  square  brackets  is  set  to  zero,  the  solution  is 


uNc  +  NT  =  I  afNc 


(C-29) 


From  Figure  C-2  it  is  clear  that  the  following  relationship  is  true 
along  the  characteristic  curves: 


dr  _  \ 

Be  "  ~N 


(C-30) 


Substituting  equation  (C-30)  into  equations  (C-28)  and  (C-29)  yields 
the  characteristic  curves  in  the  constant  n  reference  plane. 


(C-31) 
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dv  ■  i|»3  dx 

dP  -  af2  dp  -  dt 

► 

along 

dc  *  udt 

dC1  »  *5  dx 

dP  ♦  pafdu  -  +  afii  *  p«f<>2)dx 

along 

dc  *  (u  + 

o 

dP  -  pafdu  ■  (y»4  ♦  a^  ijjj  -  pa^2)dT 

along 

dc  ■  (u- 

compatibility  equations  characteristic 


af)dt 

af)dr 

curves 


Figure  C-2.  Constant  n  reference  plane  characteristic  relations 


C-32 


Equation  (C-31)  corresponds  to  the  equation  (C-28).  It  is  the  projec¬ 
tion  of  the  two-dimensional  flow  pathline  onto  the  constant  n  plane. 
Equation  (C-32)  corresponds  to  equation  (C-29).  It  is  the  projection 
of  the  two-dimensional  Mach  cone  onto  the  constant  n  plane. 

Compatibility  Equations 

The  requirement  that  the  determinant  of  the  coefficient  matrix  in 
equation  (C-26)  be  identically  zero  has  led  to  the  determination  of 
the  characteristic  curves.  Now,  substituting  those  solutions,  equa¬ 
tions  (C-28)  and  (C-29),  into  equation  (C-26)  and  solving  for  the 

J 

( j=l,2, . . . ,5)  yields  the  compatibility  equations  which  correspond  to 
the  characteristic  directions.  Specifically,  substituting  equation 
(C-28)  into  equation  (C-26)  yields  the  matrix  equation 


0  ML 


0  ML 


o  ml 


o  Ml 


o  ml 


(C-33) 


Here  the  order  of  the  coefficient  matrix  is  five  and  its  rank  is  two, 
so  there  are  three  independent  solutions.  Three  possible  solutions  are 

=  l2  =  *4  =  *5  =  0 
*  h  ■  *3  =  ‘5  =  0 
=  *2  =  *3  =  *4  *  0 


Substituting  equations  (C-34)  to  (C-36)  into  equation  (C-13)  in  turn 
yields,  respectively 


+  uvc  =  *3 

(C-37) 

PT  +  -  af2(pT  +  up^) 

*3* 

II 

(C-38) 

C1t  +  u  Cic  =  *5 

(C-39) 

dv  =  dt 

(C-40) 

dP  -  af2dP  =  ^4dt  ► 

along  d c  =  udf 

(C-41) 

dC^  =  i^gdx 

(C-42) 

Now  consider  the  characteristic  curve  given  by  equation  (C-29). 
Substituting  equation  (C-29)  into  equation  (C-26)  yields  the  matrix 
equation 


0 

0 

"A 

0 

*1 

PN? 

±afNc 

0 

0 

0 

*2 

0 

0 

±af\ 

0 

0 

l3 

0 

N?/p 

0 

±aA 

0 

h 

0 

0 

0 

0 

±afNC 

h 

=  0 


(C— 43) 


The  order  of  this  coefficient  matrix  is  five  and  its  rank  is  four  so 
there  is  one  independent  solution.  Equation  (C-43)  yields 


(C-44) 


i 
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4 


*4 


If  £4  =  1,  one  possible  solution  is 

2 

£3  =  £5  =  0;  Jlj  =  ;  Jl2  *  +  paf;  £4  =  1  (C-45) 

Substituting  equation  (C-45)  into  equation  (C-13)  yields 

2 

af  (pt  +  up^  +  pu^  -  ^j)  ±  paf(uT  +  uu?  +  P^/p  -  ip2) 

(C-46) 

+  Pt  +  uP^  -  af2(px  +  up^)  -  ip4  =  0 
or 

2 

dP  +  pa^du  =  Op 4  +  a^  +  pa^Jdx  along  d£=(u+a^)dx  (C-47) 

2 

dP  -  pafdu  =  OP4  +  ipj  -  paftp2)dx  along  d*;=(u-af)dx  (C-48) 

3.  CONSTANT  5  REFERENCE  PLANE 

For  the  constant  e  reference-plane  scheme,  all  derivatives  with 
respect  to  c  in  the  governing  equations  are  approximated  and  treated  as 
forcing  terms.  The  development  which  follows  is  entirely  analogous  to 
the  derivation  of  the  constant  n  reference-plane  relations  in  the  pre¬ 
vious  section  of  this  appendix.  Therefore,  the  corresponding  steps  in 
the  development  of  the  equations  will  all  be  shown,  but  the  arguments 
will  be  abbreviated.  Moving  all  derivatives  with  respect  to  q  to  the 
right-hand  side  of  the  governing  equations  yields 


+  vp^  +  pau^  +  p$v^  =  1|>J 

(C-49) 

+  7un  +  aP^/P  =  <P2 

(C-50) 

+  7v0  +  6Pn/p  =  *3 

(C-51) 

4 
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Px  +  vPn  ‘  af  (pt  +  Vpn}  =  *4 


(C-52) 


C1t  +  vCin  =  ^5 


(C-53) 


where 


^  =  -up^.  -  pu^  -  epv/n 


*2  ‘  'uuc  -  Vp 


*3  =  ~uv^ 


*4  =  -UP,,  ♦  af  upc+  *k 


*5  *  "uCii:  +  Vp 


Characteristic  Curves 


Forming  a  linear  combination  of  equations  (C-49)  to  (C-53)  yields 
t1(pT  +  vpn+  potun  +  pBv^  -  +  l2(uT  +  vun+  aP^/p  -  \\>2) 

+t3(vT  +  Vvn  +  BP^/p  -  ip3)  +  £4[Pt  +  7Pn  -  af2(pT  +  7pn)  -  i|^4] 


+«-5(Cix  +  vC.n  -  l>5)  =  0 


(C-54) 


Rearranging  equation  (C-54)  gives 


pn(v*l  "  af  ^4)  +  PT^j  "  af  ^4^ 


+  u  (pa^,  +  wH0)  +  u_(&,) 
n  i  t  <■  c 

+  vn(pB  Aj  +  v*3)  +  vt(*3) 


(C-55) 


+  Pn(a^2/p  +  e*3/p  +  v^4)  +  Pt(£4) 

+  Cin  ^v^5^  +  =  Vl  +  *2*2  +  £3*3  +  *4*4  +  *5*5 
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Now,  determine  if  the  £.  (j-1,2, . . . ,5)  can  be  chosen  so  that  the  W- 

J  j 

are  all  linearly  independent.  If  is  the  characteristic  normal,  then 


N-W.  =  0  (j=l,2,...,5).  This  gives 

J 

(v&j  -  a^  v&4 ) -  a^  =  0 

(C-62) 

(patj  +  vt2)Nri  +  t2NT  =  0 

(C-63) 

(pBt1  +  v£3)Nn  +  £3Nt  =  0 

(C-64) 

(a£2/p  +  8*3/p  +  vt4)Nri  +  £4Nt  =  0 

(C-65) 

(Vt5)Nn  +  t5NT  =  0 

(C-66) 

or,  in  matrix  form: 


’'Vt 

0 

0 

-af(vNn+Nx) 

0 

pcNn 

vN+N 

n  x 

0 

0 

0 

p8N 

n 

0 

*VNx 

0 

0 

0 

aNn/p 

8Nn/p 

vN  +N 

n  x 

0 

0 

0 

0 

0  vN 

+N 

(C-67) 


Setting  the  determinant  of  the  coefficient  matrix  to  zero  yields 

(7Nn  +  NT)3[(7Nn  +  Nt)2  -  af2Nn2(a2  +  $2)]  =  0  (C-68) 

Equation  (C-68)  has  two  possible  solutions: 


v  Nn  +  Nt  =  0 


v  N  +  N  =  ±  a^N  a* 
n  x  f  n 


(C-69) 
(C— 70) 


where 


a*  =  (a2  +  (32)*. 


Noting  from  Figure  C-3  that  dx/dri  =  -N^/N^.,  equations  (C-69)  and  (C-70) 
can  be  rewritten  as 


dT  V 


&  -  v  ?  af 


(C-71) 

(C-72) 


respectively.  Equations  (C-71)  and  (C-72)  are  the  characteristic  curves 
in  the  constant  x,  reference  plane. 


I 
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3du  -  adv  *  ~ 

2 

dP  -  dp  =  <J^dx 
dC^  *  ij/gdx 


paaf  p3af 

dP  +  — J-  du  +  t~  dv 
a* 


a’ 


?  paaf  pBa-ij;- 
(*4  +  af  +  -&r  $2  +  r -)<Jt 


paa-  p3af 

dP - _L  du  +  <jv  * 

a*  a* 


(*4  +  «f  4»j  "  — 


paaf^2  pBafi|/3 


a 


^ — —)  dx 


along  dn  =  vdx 


along  dn  *  (v  +  a^Jdx 


along  dn  s  (v  -  a^a*)dx 


compatibility  equations  characteristic  curves 


Figure  C-3.  Constant  c  reference  plane  characteristic  relations 


Consider  the  characteristic  curve  given  by  equation  (C-69). 
stituting  equation  (C-69)  into  equation  (C-67)  yields 


0 

0 

0 

0 

0 

*1 

p“Nn 

0 

0 

0 

0 

£2 

pSNn 

0 

0 

0 

0 

Z3 

0 

aNn/p 

3N  /p 
n  K 

0 

0 

U 

0 

0 

0 

0 

0 

A5 

(C-73) 


Three  possible  solutions  to  equation  (C-73)  are 

£  2  =  1  >  ^3  =  -a/3 

-  £ 2  -  £3  =  £g  =  0;  £4  =  1 
£l  =  £2  =  £3  =  £4  =  °;  £g=l 


(C-75 

(C-76 


Substituting  equations  (C-74)  to  (C-76)  in  turn  into  equation  (C-55) 
yields,  respectively 

v  ^  +  u^  +  aP^/p  +  v(-a/6)vr|  -  a/3  +  3/p(-a/3)Pn 

(C-77 


Vv  V  af2pT  +  vPn  +  p,  -  *4 
V  Cin  +  c1t  =  *5 


=  <l»2  " 


(C-78 

(C-79 


which  may  be  written  as 


(C-8( 


gdu  -  adv  =  (6^2  “ 

dP  -  a^dp  =  tp^di 

dC.  =  ipc  dx 
1  □ 


along  dn=vdx 


(C-81 


( C-82 


Now  consider  the  characteristic  curve  given  by  equation  (C-70). 
Substitution  of  this  result  into  equation  (C-67)  yields 


±afNna* 

0 

0  -a 

f<±afNn“* 

)  0 

*1 

paNri 

±a,N  a* 
f  n 

0 

0 

0 

h 

P3Nn 

0 

±afNna* 

0 

0 

l3 

0 

aNn/P 

BNn/p 

±afNna* 

0 

*4 

0 

0 

0 

0 

±a^N  a* 
f  n 

*5 

Equation  (C-83)  has  one  independent  solution: 

p 

h  =  af  V  h  =  B  *3’  aiZ  +  eA3  =  +  Pafa%;  %  =  0  (C-84 

If  =  1  is  chosen,  then 
2 

*4  =  af  *2  =  +  paa^/a*;  JI3  =  +  p^a^/a*; 

(C-85 

A4  =  1 »  %  =  0 

Substituting  equation  (C-85)  into  equation  (C-55)  yields  the  compati 
bility  equations 

paaf  psaf  2  P“af  P^f 

dP  +  — du  +  —5-  dv  =  (i|)4  +  af  ^  +  —jt  ip2  +  — ^3)dT 

( 

along  dn  =  (v  +a^x*)dT 
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APPENDIX  D 


THE  SPECIES  CONTINUITY  EQUATION  INTEGRATION  SCHEME 

1.  STIFF  EQUATIONS  AND  IMPLICIT  VS.  EXPLICIT  INTEGRATION  METHODS 

It  is  well  known  that  integration  of  the  species  continuity  equa 
tions  for  nonequilibrium,  chemically  reacting  flow,  when  the  flow  is 
near  equilibrium,  requires  special  treatment  because  the  equations 
become  quite  "stiff"  (Refs.  3,  5,  and  31).  In  order  to  understand 
the  concept  of  a  stiff  differential  equation,  consider  the  following 
general  form  for  an  ordinary  differential  equation. 

%  =  f(x,y)  (D-l) 

The  existence  of  a  unique  solution  to  this  equation  requires  that 
f(x,y)  must  be  continuous  and  that  it  satisfy  a  criterion  known  as  the 
Lipschitz  condition  (Ref.  40).  The  Lipschitz  condition  is  stated  in 
terms  of  the  Lipschitz  constant,  L,  which  is  defined  by 

l  =  |  |  (D-2) 

oj 

If  the  ordinary  differential  equation  is  stiff,  it  will  have  a  large 
Lipschitz  constant  while  its  solution  behaves  like  a  polynomial  (i.e., 
the  solution  has  little  exponential  growth  or  decay). 

The  relevance  of  the  Lipschitz  constant  to  the  species  contin¬ 
uity  equation  can  be  understood  from  the  following  discussion.  The 


species  continuity  equation  can  be  written  in  functional  form  as 


DC .  0. 

=  —  =  f^P.T.C^ )  (i=l,...,n)  ( D-3) 

Note  that  f^ (p,T,C^ )  in  equation  (D-3)  has  dimensions  of  divided  by 
time  and  therefore,  the  partial  derivative  (3-f ^/3C^ )p  j  must  have 
dimensions  of  (time)"^.  Following  Ref.  (25),  a  local  characteristic 
time  for  the  rate  process  can  be  defined  by 


/s 

T 


1 

(9^/30.) 


P  »T 


(D-4) 


Now,  equation  (D-3)  can  be  written  as 


DC.j  x  (p  »T  ,C. ) 
"DT  =  x(P  ,T,C. 


(i«l,  —  ,n) 


where 


X(P,T,C1) 


f< 

<3V5Ci'p,t 


(D-S) 


(D-6) 


Using  equation  (D-3)  and  the  definition  of  the  Lipschitz  constant 
(D-2),  the  Lipschitz  constant  for  the  species  continuity  equation 
becomes 


L  = 


3f i (P»T»Ci )p,T 
- 5C. -  1 


(D— 7) 


Consider  the  change  in  according  to  equation  (D-3)  for  an  increment 
of  time  equal  to  the  local  relaxation  time. 


DCi 

T5T 


C2"C1 


(D-8) 
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In  equilibrium  flow,  the  local  characteristic  time,  x,  for  the  chemi- 

DC 

cal  reaction  approaches  zero  and  since  must  remain  finite  for  any 
real  physical  process,  the  change  in  mass  fraction  also  approaches 
zero.  Thus,  the  Lipschitz  constant  becomes  large  while  the  solution 
behaves  like  a  polynomial  and  the  differential  equation  is  stiff.  For 
the  case  of  frozen  flow,  the  mass  fraction  does  not  change  and  x 
approaches  infinity,  so  DC./Dt  goes  to  zero. 

Many  numerical  techniques  have  been  proposed  for  the  solution  of 
the  species  continuity  equations  in  near  equilibrium  flow.  Among  these 
techniques  are  linearized  approaches,  conversion  of  the  differential 
equations  to  integral  equations,  modified  Range-Kutta  techniques, 
higher-order  predictor-corrector  methods,  and  implicit  methods.  Curtiss 
and  Hirshfelder  (32)  and  Seinfield  (33)  provide  useful  reviews  of  the 
various  numerical  methods  that  have  been  used  for  stiff  differential 
equations.  Cline  (41)  analyzed  and  tested  a  number  of  proposed  schemes 
in  his  analysis  of  three-dimensional,  steady,  nonequilibrium  flow. 

It  is  possible  to  categorize  the  numerical  methods  he  tested  into  two 
groups:  explicit  and  implicit  methods.  In  an  explicit  method,  the 
solution  at  the  unknown  point  is  expressed  entirely  in  terms  of  known 
points,  while  in  an  implicit  method,  the  solution  at  the  unknown  point 
is  expressed  in  terms  of  both  the  known  points  and  the  unknown  point. 
Cline  concluded  that  explicit  schemes,  and  predictor-corrector  methods, 
were  stable  only  for  very  small  step  sizes  (which  would  lead  to  exces¬ 
sive  computational  time).  Tyson  and  Kliegel  (42)  show  that  a  step 
size  of  the  order  of  the  relaxation  time  x  is  necessary  to  insure 
stability  when  the  first-order  Euler  method  (explicit)  is  applied  to 
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the  species  continuity  equations.  Cline  also  found  that  only  sane  of 
the  implicit  methods  he  tested  gave  adequate  results.  He  concluded 
that  a  method  proposed  by  Lomax  and  Bailey  (34)  provided  the  best 
results.  It  is  this  technique,  which  could  be  called  a  second-order, 
implicit  Taylor  expansion,  which  is  used  in  this  analysis. 

2.  DERIVATION  OF  THE  SECOND-ORDER,  IMPLICIT,  TAYLOR  EXPANSION  SCHEME 
To  illustrate  the  derivation,  consider  the  ordinary  differential 
equation: 


^  =  f(x,y)  (D-9) 

Here,  f(x,y)  may  be  nonlinear  in  y  and  y(x)  is  assumed  to  be  continuous 
and  differentiable  so  it  can  be  expanded  in  a  Taylor  series.  Expanding 
through  second-order  terms  yields 


yn+l  =  *n  +  dx !n  (xn+l 


v  )  +  I  (xn+l"xn^2 

V  dx2  >n  2 


+  0Ux3) 


(D-10) 


where 


ti  ■*&... 


and  n  denotes  the  nth  mesh  point.  Note  that 


0 1  =  &  (3x! }  =  &  Cfn(x*y)] 

(be  n  n 


(D-U) 
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(D- 12) 


3f  3f  ,u  3f_ 

_ n  +  _ n  dy  _  _ n  + 

3x  3y  dx  ~  3x 


!!a 

dy 


Substituting  equation  (D-12)  into  equation  (D-10)  yields 
»■»!  *  »„  +  hfn  *  T  hiT  +  f„  V1  +  0<h3) 


(0-13) 


where  h  =  xp+1  -  xn-  Now,  fn  inside  the  square  brackets  in  equation 
(D-13)  is  replaced  as  follows 


f 


n 


^n+1  ~  _  ^n+1  ~  ^n 

xn+l  '  xn  '  h 


(0-14) 


Substituting  equation  (D-14)  into  equation  (D-13)  yields  the  desired 
result. 


2  3f_ 


af„  (y 


Vl  ‘  +  hfn  *  T  C1T  +  ly1 


n+1 


-  yn)  o 

: - — ]  +  0  ( h3 )  ( D— 15) 


The  derivation  for  the  species  continuity  equations  is  entirely 
analogous  to  the  above  derivation.  The  species  continuity  equations 
are  restated  here  for  convenience. 

DC.  o. 

Ut1  = =  f^P.T.C.)  (1=1 . ft)  (D-16) 

Note  that  equation  (D-16)  is  a  first-order,  ordinary  differential  equa¬ 
tion  when  applied  along  a  particle  path  in  unsteady  flow  or  along  a 
streamline  in  steady  flow.  Along  the  one-dimensional  trajectory,  the 
substantial  derivative  can  be  written  as  D(  )/Dt  =  Vd(  )/ds  where  s 
is  an  elemental  path  length  and  V  Is  the  velocity  magnitude.  Thus, 
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equation  (D-16)  becomes 


dC.  o- 

-fe-  =  -J  =  gi(p,v»T,ci)  0=1,. ...ft) 


(0-17) 


The  i  subscript  is  temporarily  dropped  in  the  following  development  for 
clarity.  Expanding  C  in  a  Taylor  series  through  second-order  terms 


yields 


0  =  c  ♦  4£|  (S  .  s  )  +  A  I  'Vr5.1' 

Ln+1  Ln  ds  I  'sn+l  V  .2  ' n 

n  ds 


(D-18) 


+  0  (As*5 ) 


—  (—1  ) 
ds  'ds  V 


£[gn(p,v,T,c)] 


( D- 19 ) 


89n  do  "n  dV  ,  "n  dT  ,  3gn  dC 
3p  ds  9V  ds  dT  ds  aC  ds 


(D-20) 


Note  that  the  last  term  in  equation  (D-20)  will  become  a  summation  over 
all  species  when  the  i  subscript  is  re-introduced.  Substituting  equa¬ 
tion  (D-20)  into  equation  (D-18)  yields 


n  dV  ±  dyn  dT  ,  ayn 


r  -  r  x  Ur,  x  r  "  u  .  n  u,  .  _ ri  ui  .  m  _  n  " 

Cn+1  "  Cn  hgn  '  a0  ds  aV  ds  aT  ds  aC  gnJ  2 


i] 


m 


1* 


where  h  =  s„^,  -  s„.  Now,  let  ds  -  sn+1  -  sp;  dT  =  Tn+1  -Tn,  etc.. 


n+1  Jn 


n+1  n! 


and  let  gn  inside  the  square  bracket  be  replaced  by 

„  _Cn+l-Cn 
% - h - 


Substituting  into  equation  (D-21)  yields 


Cn+1  =  Cn  +  hgn  +  2  ^  3p  ^pn+l  “  pn^  +  3V  ^Vn+1  "  V 


39, 


39 


39 


Re-introducing  the  i  subscript  yields 


Kin 


.  h 
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3gin 
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where 


K.  =  Ci  . ,  =  Ci  . 
in  n+l  n 


(D-22) 


+  if  (Tnn  -  Tn>  +  if  <Cn+l  ‘  Cn»  *  °<h  >  O’23) 


(D-24) 


The  partial  derivatives  in  equation  (D-24)  are  determined  analytically 
When  expanded,  equation  (D-24)  becomes  a  system  of  h  simultaneous, 
linear,  algebraic  equations  which  is  easily  solved  by  standard  methods 
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3.  ERROR  CONTROL 

The  implicit  scheme  derived  in  the  previous  section  has  been  found 
to  be  stable  for  even  large  step  sizes  in  both  Cline's  work  (41)  and  in 
the  subject  research.  However,  the  method  is  subject  to  truncation 
error  and  thus  consideration  of  step  size  control  must  be  given. 

Equation  (D-24)  is  of  the  form 

^  +  0(h4)  1D-25) 

ds  n 


"computed  term"  "third  order 

term" 


The  objective  in  the  error  control  scheme  is  to  estimate  the  third- 
order  term  and  require  that  the  ratio  of  this  term  to  the  computed 
term  be  less  than  a  specified  tolerance.  The  third-order  term  is  esti 
mated  as  follows. 
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(D-26) 


then 
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and  RATIO  is  computed  for  each  of  the  chemical  species.  If,  after  the 
last  intermediate  integration  step,  the  maximum  /a lue  of  RATIO  is  less 
than  a  user  specified  tolerance  (TOL),  then  the  integration  proceeds. 

If  the  required  tolerance  is  not  achieved,  then  the  number  of  inter¬ 
mediate  points  is  doubled  and  the  integration  procedure  is  restarted 
from  L.  Also,  if  the  maximum  value  of  RATIO  is  two  orders  of  magni¬ 
tude  less  than  TOL,  then  NINT  is  -educed  by  one-half  before  the  inte¬ 
gration  from  L+l  to  L+2  proceeds.  Figures  D-2  and  D-3  show  the  effect 
of  NINT  in  integrating  between  given  end  conditions  for  equilibrium 
and  nonequilibrium  conditions  at  L.  In  the  first  two  sections  of  the 
tables  in  those  figures,  NINT  has  been  specified,  while  in  the  last 
section,  TOL  has  been  specified  and  NINT  has  been  determined  by  the 
error  control  procedure  described  above.  Note  that  in  both  figures, 
failure  to  use  intermediate  points  yields  poor  results  for  the  species 
mass  fractions  and  However,  for  the  cases  shown  here,  even  a  small 
number  of  intermediate  points  improves  the  results  significantly. 
Computational  time  for  the  subject  problem  is  almost  directly  propor¬ 
tional  to  NINT.  Experience  indicates  that  a  loose  tolerance  (5  to  10%) 
should  be  used  until  the  solution  is  near  convergence  and  then  TOL  can 
be  reduced  to  the  desired  degree. 

Figure  D-4  illustrates  the  effect  of  variable  property  gradients 
between  L  and  L+l  on  the  number  of  intermediate  points  required  to 
achieve  a  tolerance  of  one  percent.  Part  (a)  of  the  figure  corres¬ 
ponds  to  equilibrium  conditions  at  L.  Case  1  indicates  that  for  small 
gradients  in  density  and  temperature,  the  integration  from  equilibrium 
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Figure  D-2.  The  effect  of  NINT  intermediate  points  -  equilibrium 
initial  conditions  (H-F  system). 
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Figure  D-3.  The  effect  of  NINT  intermediate  points 
initial  conditions  (H-F  system). 
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conditions  may  require  a  large  number  of  intermediate  points  to  achieve 
a  one  percent  tolerance.  However,  cases  2  and  4  show  that  even  modest 
gradients  in  temperature  in  density  reduce  NINT  dramatically.  In  com¬ 
paring  cases  2  and  3  note  that  NINT  is  not  sensitive  to  the  sign  of 
the  temperature  gradient.  Case  5  corresponds  to  tripling  the  relaxa¬ 
tion  time  between  L  and  L+l. 

In  part  (b)  of  figure  D-4,  the  conditions  at  L  are  not  in  equili¬ 
brium.  The  gradients  represented  by  case  1  are  the  same  as  those  that 
existed  between  L-l  and  L  in  the  data  from  which  this  information  was 
taken.  Cases  2  and  3  show  that  NINT  is  sensitive  to  a  change  in  the 
temperature  gradient  but  that  it  is  not  affected  by  a  change  in  the 
density  gradient.  Cases  4  and  5  illustrate  the  effect  of  the  velocity 
gradient  (or  relaxation  time)  on  NINT.  In  general,  discontinuous 
gradients  in  temperature  and  velocity  between  the  L-l  to  L  interval 
and  the  L  to  L+l  interval  will  require  a  higher  value  of  NINT  to  achieve 
the  specified  tolerance  than  for  the  case  of  continuous  gradients. 

NINT  is  also  affected  by  the  density  gradient  but  to  a  lesser  degree 
than  for  temperature  and  velocity.  The  significant  results  concerning 
the  use  of  error  control  are  summarized  as  follows: 

1.  Most  of  the  computational  time  required  to  solve  the  subject  problem 
is  related  to  the  integration  of  the  species  continuity  equations. 
Therefore,  the  total  execution  time  is  nearly  directly  proportional 
to  NINT. 

2.  Failure  to  use  intermediate  points  in  integrating  the  species  con¬ 
tinuity  equations  will,  in  general,  yield  poor  results  for  species 
mass  fractions  and 
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3.  Integrating  from  equilibrium  conditions  with  very  small  gradients 
in  temperature  and  density  and  a  small  value  of  TOL  can  yield  very 
large  numbers  of  intermediate  points. 

4.  For  equilibrium  initial  conditions,  NINT  is  not  particularly  sen¬ 
sitive  to  the  size  or  magnitude  of  temperature,  density,  or  velo¬ 
city  gradients  as  long  as  the  conditions  in  conclusion  3  do  not 
occur. 

5.  For  nonequilibrium  initial  conditions,  NINT  is  sensitive  to  grad¬ 
ients  in  temperature  and  velocity  while  it  is  only  mildly  sensitive 
to  the  density  gradient.  If  the  property  gradients  between  L  and 
L+l  do  not  match  corresponding  gradients  between  L-l  and  L,  then 
more  intermediate  points  are  likely  to  be  required. 

6.  Conclusions  1  to  5  appear  to  be  valid  for  the  several  chemistry 
systems  investigated  in  this  research. 
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TREATMENT  OF  SUBSONIC  INLET  POINTS 

1.  IMPORTANCE  OF  COMPUTATIONAL  BOUNDARY  CONDITIONS 

Roache  (43)  remarks  that  all  the  flow  patterns  of  common  gases 
and  liquids  are  solutions  to  the  same  partial  differential  equations. 
The  f 1 ows  are  distinguished  only  by  boundary  and  initial  conditions 
and  by  iimilarity  parameters  like  the  Reynolds  number.  He  concludes: 
"It  is  therefore  not  surprising  that  the  specification  of  computa¬ 
tional  boundary  conditions,  besides  affecting  numerical  stability, 
greatly  affects  the  accuracy  of  the  FDE  solution."  Cline  (23)  attri¬ 
butes  long  computational  times  in  time-dependent  solutions  of  con¬ 
verging-diverging  nozzle  problems  to  poor  treatment  of  the  boundaries 
and  inefficient  schemes  for  the  interior  points.  Moretti  (44)  and 
Abbett  (45)  have  shown  that  reflection,  extrapolation,  and  one-sided 
difference  techniques  give  poor  results  when  applied  to  solid  wall 
boundaries.  Cline  (4,  23)  has  used  a  reference-plane  characteristic 
method  for  boundary  points  in  time-dependent  analyses  of  nozzle 
problems.  This  technique  is  also  used  at  all  boundary  points  in  the 
subject  research. 

2.  CHOICE  OF  INLET  BOUNDARY  CONDITIONS 

At  the  inlet,  a  constant  n  reference-plane  characteristic  method 
is  used  and,  for  subsonic  flow,  only  one  characteristic  curve  lies 


within  the  flow  field  (see  Appendix  C).  Since  there  are  4+n  unknowns 
at  the  inlet  (u,  v,  P,  p,  and  ,  i=l,...,n)  and  only  one  character¬ 
istic  relation,  3+n  additional  conditions  must  be  specified. 

The  "correct"  specification  of  steady  flow  inlet  boundary  con¬ 
ditions  is  a  difficult  problem;  a  generally  preferred  set  of  condi¬ 
tions  is  not  known  at  the  present  time.  It  is  expected  that  the 
combustor  analysis  will  be  performed  independently  of  the  nozzle 
analysis.  The  combustor  analysis  may  or  may  not  include  two- 
dimensional  effects.  In  any  case,  the  distributions  of  the  fluid 
dynamic  properties  (u,  v,  p,  and  p)  at  the  combustor  exit  may  not 
provide  an  appropriate  set  of  boundary  conditions  for  the  given  nozzle 
problem.  The  total  conditions  at  the  combustor  exit,  however,  are 
global  properties  of  the  flow  which  must  be  conserved  between  the 
combustor  exit  and  the  nozzle  inlet.  The  specification  of  total 
conditions  at  the  inlet  places  a  constraint  on  the  energy  of  the  flow 
but  still  allows  flexibility  in  the  solution  for  the  static  proper¬ 
ties  at  this  point. 

The  specification  of  total  enthalpy  HQ  is  appropriate  for  the 
subject  problem.  Equation  (A-3)  shows  that  the  total  enthalpy 

p 

(h  +  V  /2)  must  be  constant  along  streamlines  in  a  steady  flow. 
Therefore,  the  total  enthalpy  distribution  at  the  inlet  is  required 
to  match  that  at  the  combustor  exit.  Cline  (4,  23)  uses  the  inlet 
flow  angle  e  as  a  boundary  condition.  It  has  also  been  used  for  the 
subject  problem.  The  flow  angle  may  be  determined  by  experiment  or 
from  the  following  procedure: 


1.  Determine  representative  values  of  the  ratio  of  frozen  specific 

heats  and  the  gas  constant  at  the  combustor  exit. 

2.  Generate  a  "long  inlet"  geometry  by  adding  six  to  ten  mesh  points 
upstream  of  the  nozzle  inlet,  simulating  a  constant-area  duct. 

3.  Using  the  program  of  the  subject  research  for  isentropic  constant 

specific  heat  ratio  flow,  compute  the  solution  to  the  long  inlet 

problem  described  in  the  preceding  two  steps.  Inlet  boundary 

conditions  for  the  long  inlet  problem  are:  a  uniform  distribution 

of  zero  flow  angle  9,  and  total  pressure  PQ  and  total  temperature 

T  based  on  combustor  exit  conditions, 
o 

4.  The  flow  angles  at  the  nozzle  inlet  for  the  converged  long  inlet 
problem  can  be  used  as  the  inlet  flow  angles  for  the  reacting 
flow  problem. 

Note  that  the  arbitrary  specification  of  zero  flow  angle  at  the  nozzle 
inlet  has  been  found  to  produce  unreasonable  velocity  profiles  at  the 
nozzle  inlet  for  some  nozzle  geometries. 

A  uniform  distribution  of  static  pressure  across  the  combustor 
exit  is  one  of  the  results  expected  from  the  combustor  analysis.  The 
specification  of  9,  HQ,  C.  and  a  uniform  static  pressure  distribution 
was  tested  as  an  inlet  boundary  condition  set.  For  the  cases  tested, 
this  produced  a  velocity  distribution  which  had  a  range  from  very  high 
values  at  the  center  to  near  zero  values  at  the  wall  after  only  20  to 
30  time  planes.  The  static  pressure  distribution  at  the  inlet  from 
the  "long  inlet"  solution  yielded  reasonable  velocity  distributions. 

The  set  of  inlet  boundary  conditions  which  has  been  used  in  the 
subject  research  is  the  specification  of  9,  HQ,  C^,  1*1,..., n  and  the 


total  pressure  PQ  based  on  combustor  exit  conditions.  Here,  PQ  is  the 
pressure  that  would  exist  as  a  result  of  an  isentropic  stagnation  from 
combustor  exit  conditions  with  frozen  composition. 

3.  RELAXATION  OF  THE  SPECIES  MASS  FRACTIONS 

In  using  total  enthalpy  Hq  and  total  pressure  PQ  as  inlet  boun¬ 
dary  conditions  in  the  constant  n  reference-plane  scheme,  it  is  found 
that  radial  distributions  of  density,  temperature,  and  velocity  develop 
which  do  not  necessarily  match  similar  distributions  at  the  combustor 
exit.  Therefore,  a  hypothetical  particle,  upon  making  the  transition 
from  the  combustor  exit  to  the  nozzle  inlet,  encounters  discontinuities 
in  the  flow  field.  If  not  properly  accounted  for,  these  discontinui¬ 
ties  provide  inconsistent  kinetics  data  at  the  nozzle  inlet  from  which 
integration  of  the  species  continuity  equations  proceeds. 

A  series  of  numerical  experiments  was  performed  to  study  the 
effects  of  the  discontinuities  described  above.  Two  different  chem¬ 
istry  systems,  in  different  nozzles,  and  in  various  degrees  of  "near 
equilibrium"  have  been  subjected  to  discontinuities  of  different  types 
and  magnitudes.  The  geometries  and  chemistry  systems  are  presented 
in  Figure  E-l. 

Consider  a  particle  at  the  combustor  exit.  Note  that  in  func- 


tional  form  (see  Appendix  A) 

o.j  -  o.j  ( o.T ,C.. )  (i-1, . . .  ,n) 

(E-l) 

=  ( T , o.  *C.j )  =  p»T 

(E-2) 
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C-H-O-N  GEOMETRY 

H20  +  M  t  OH  +  H  +  M 

H2  +  M  t  2H  +  M 

02  +  M  ?  20  +  M 

C02  +  H  t  OH  +  CO 

H20  +  H  t  H2  +  OH 

H20  +  0  4  20H 

OH  +  H  t  H2  +  0 

OH  +  0  *  02  +  H 

C-H-O-N  REACTION  MECHANISM 

Figure  E-l. 


H-F  GEOMETRY 

F2  +  M  t  2F  +  M 

HF  +  M^H  +  F  +  M 

H2  +  M  t  2H  +  M 

HF  +  F  J  F2  +  H 

HF  +  H  t  H2  +  F 

2HF  t  F2  +  H2 

H-F  REACTION  MECHANISM 


C-H-O-N  and  H-F  systems. 


Since  all  the  static  properties  and  the  species  mass  fractions  are 
known  at  the  location  of  the  particle,  the  species  source  functions 
^  and  the  energy  source  term  ^  can  be  computed  at  this  point.  In  a 
real  flow,  there  can  be  no  discontinuities  of  the  static  properties 
and  species  mass  fractions  at  the  junction  between  the  combustor  exit 
and  the  nozzle  inlet.  Therefore,  from  equations  (E-l)  and  (E-2),  the 
distributions  of  o.  and  ^  along  the  flow  direction  must  also  be  con¬ 
tinuous  at  this  junction.  The  discontinuities  do  not  represent 
reality.  They  are  merely  the  consequence  of  an  attempt  to  link 
together  separate  combustor  and  nozzle  analyses.  Their  effect  should 
not  be  felt  within  the  flow,  but  should  be  manifest  at  the  inlet  mesh 
points  before  the  species  continuity  equations  are  integrated  through 
the  flow  field.  This  is  accomplished  by  allowing  the  species  mass 
fractions  at  the  combustor  exit  to  adjust  to  the  static  conditions  at 
the  nozzle  inlet.  The  discontinuities  are  replaced  by  ramp  functions 
outside  (just  before)  the  inlet  mesh  points  as  shown  in  Figure  E-2. 

The  discontinuities  in  density  and  temperature  become  continuous, 
linear  transitions  from  combustor  to  inlet  conditions  over  the  period 
of  a  variable  "relaxation  time."  The  time  is  chosen  so  that  the  change 
in  the  energy  source  termip^  is  minimized  between  the  combustor  exit 
and  the  nozzle  inlet. 

Figure  E-3  shows  curves  of  ^  versus  relaxation  time  for  the  H-F 
system  starting  from  near-equilibrium  conditions.  The  effects  of 
temperature  and  density  ramps  are  shown  separately.  A  zero  relaxation 
time  corresponds  to  a  discontinuity  at  the  junction  between  the 


Figure  E-2.  Unsteady  relaxation  at  the  inlet 


initial  conditions  (H-F  system). 


combustor  exit  and  the  nozzle  inlet.  Note  that  a  temperature  discon¬ 
tinuity  of  just  -20°R  changes  the  value  of  ^  by  several  orders  of 
magnitude  (from  2x  10^  to  974 x  10^).  The  size  of  the  temperature  and 
density  discontinuities  in  Figure  E-3  are  typical  for  the  cases 
studied  during  this  research.  Also,  note  that  with  increasing  relax¬ 
ation  time,  the  values  of  ^  approach  zero.  This  is  because  the  chem¬ 
ical  systems  move  toward  the  equilibrium  condition  where  the  species 
source  functions  (and  consequently  the  energy  source  term)  are  zero. 
Table  E-I  shows  the  effect  of  relaxation  time  on  the  species  mass 
fractions.  Note  that  even  relatively  long  relaxation  times  produce 
small  changes  in  the  species  mass  fractions.  These  changes  are 
typically  two  percent  or  less  with  the  exception  of  the  change  for  F2 
which  is  present  in  minute  quantities  in  the  reactive  mixture.  Figure 
E-4  shows  versus  relaxation  time  curves  with  initial  conditions 
corresponding  to  various  locations  in  the  ODK  [Reference  (10)]  solu¬ 
tion  of  the  H-F  problem.  Note  that  for  all  sets  of  initial  conditions 
(i.e.,  various  degrees  of  near-equilibrium),  <1^  approaches  zero  for 
large  relaxation  time.  Figure  E-5  shows  ^  versus  relaxation  time  for 
the  C-H-O-N  system.  Initial  conditions  are  near  equilibrium.  Table 
E-I I  shows  the  changes  in  the  reactive  species  mass  fractions  as  a 
function  of  relaxation  time. 

There  are  several  important  conclusions  regarding  the  data  des¬ 
cribed  above: 

1.  All  curves  of  4^  versus  relaxation  time  approach  zero  as  relaxa¬ 
tion  time  increases.  This  is  because  the  chemical  systems  move 
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versus  relaxation  time  -  various  nonequilibrium  initial  conditions 
(H-F  system). 


versus  relaxation  time  -  equilibrium  initial  conditons  (C-H-O-N  system). 


toward  an  equilibrium  condition  (where  the  species  source  func¬ 
tions  go  to  zero)  with  increasing  relaxation  time. 

2.  At  high  temperatures  and  near  equilibrium  conditions,  is  quite 
sensitive  to  variations  in  temperature  and  density.  For  example, 
a  20°R  temperature  perturbation,  with  no  relaxation,  can  change 
the  computed  value  of  by  several  orders  of  magnitude. 

3.  For  the  chemistry  systems  and  static  property  gradients  studied 
in  this  research,  relatively  small  changes  in  composition  occur 
during  the  relaxation  from  combustor  exit  to  nozzle  inlet  condi¬ 
tions.  Typically,  species  mass  fractions  change  by  less  than  two 
percent. 

Conclusion  3  is  significant  with  respect  to  the  use  of  PQ  based  on 
combustor  exit  conditions  for  the  calculation  of  inlet  points.  Since 
the  composition  changes  with  relaxation  between  the  combustor  exit  and 
nozzle  inlet,  PQ  based  on  combustor  exit  conditions  is  only  an  approx¬ 
imation  at  the  inlet.  However,  due  to  the  very  small  changes  in  the 
species  mass  fractions  between  those  two  points,  it  is  a  reasonable 
approximation. 

Failure  to  relax  the  species  mass  fractions  at  the  inlet  leads 
to  numerical  difficulties  for  two  primary  reasons.  First,  without 
relaxation,  the  data  set  (p,T,C..,  1=1,..., n)  at  the  inlet  is  inconsis¬ 

tent.  As  noted  above,  discontinuities  can  cause  the  energy  source  term 
^  to  differ  by  several  orders  of  magnitude  from  its  value  at  the  com¬ 
bustor  exit.  This  causes  a  "hard  start"  condition  in  the  integration 
of  the  species  continuity  equations  away  from  the  inlet.  Second,  if 


the  species  continuity  equations  are  not  relaxed  before  the  inlet, 
they  will  relax  within  the  computed  flow  field.  This  can  produce 
distributions  of  4^  within  the  flow  with  dramatic  oscillations  near 
the  inlet.  This  4^  distribution,  in  turn,  distorts  the  flow.  Figure 
E-6  illustrates  the  difficulties  described  above.  The  "total  enthalpy 
error"  in  that  figure  is 


H  Error 
o 


’  H0R 


2 


(E-3) 


where  HqL  is  the  total  enthalpy  at  station  L,  HqR  is  the  total  enthalpy 
at  the  combustor  exit,  and  V£  and  Vj  are  the  velocity  magnitudes  at 
the  nozzle  exit  and  inlet  respectively.  Shown  in  Figure  E-6  is  a  plot 
of  total  enthalpy  error  and  4>k  at  each  centerline  L  station  in  the  H-F 
problem  for  various  relaxation  times  at  the  inlet.  Note  that  for  very 
short  relaxation  times,  4>k  does  not  have  time  to  adjust  to  conditions 
within  the  flow  and,  as  a  consequence,  the  rapid  adjustment  of  4^  at 
the  inlet  can  produce  large  total  enthalpy  errors.  Once  this  error 
is  established  at  the  inlet,  no  recovery  is  made. 

The  proper  relaxation  time  is  determined  as  part  of  the  overall 
numerical  algorithm.  The  scheme  is  illustrated  in  Figure  E-7  which 
shows  a  typical  4>k  versus  relaxation  time  curve.  In  that  figure,  4JkR 
is  computed  for  the  conditions  at  the  combustor  exit,  4>kQ  is  computed 
for  the  species  mass  fractions  at  the  combustor  exit  and  density  and 
temperature  at  the  nozzle  inlet  with  zero  relaxation  time,  and 
is  computed  for  the  same  conditions  as  4>kQ  but  with  a  relaxation  time 
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enthalpy  error  and  profiles  for  various  relaxation  times 


At.  The  values  of  if/^R  corresponding  to  each  inlet  mesh  point  are  com¬ 
puted  as  part  of  the  initialization  procedure  for  the  overall  algorithm. 
They  are  only  computed  once  and  then  are  stored.  Consider  time  level 
N  in  the  solution  of  the  subject  problem,  is  computed  at  each 
streamline  origin  for  the  static  conditions  prevailing  at  the  inlet. 
There  is  a  minimum  relaxation  time  and  if  |i^q|  <  |g^R|  ,  then  the 
minimum  relaxation  time  is  used.  This  is  because  is  known  to 
approach  zero  with  increasing  relaxation  time.  Otherwise,  the  species 
continuity  equations  are  relaxed  through  time  At^  (see  Figure  E-7) 
which  is  the  relaxation  time  computed  for  a  given  streamline  during 
previous  time  steps.  Then  the  species  continuity  equations  are  inte¬ 
grated  along  the  given  streamline  to  the  nozzle  exit.  Next,  a  check 
is  made: 


K&t-'tW  I 

'*k0  *  *kRl 


.05 


(49) 


If  the  check  is  satisfied,  then  A^  is  not  changed  and  it  is  saved  for 
use  during  the  next  time  step  in  the  overall  algorithm.  However,  if 
the  check  is  not  satisfied,  then  a  straight  line  is  extended  from  <J^q 
through  to  intersect  i|^R.  The  relaxation  time  Atg  corresponding 

to  that  intersection  is  used  during  the  next  time  step.  The  same 


procedure  continues  for  subsequent  time  steps  as  shown  in  Figure  E-7 
until  the  check  is  satisfied. 
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